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1.  Abstract 


HyPerComp  is  developing  advanced  computational  tools  to  model  high  speed  flows  with  MHD  effects. 
Unstructured,  multi-dimensional  hypersonic  MHD  codes  have  been  developed  at  HyPerComp  to  study 
supersonic  viscous  flows  in  a  self-consistent,  fully  coupled  manner.  Effects  of  thermal,  chemical  and 
internal  mode  dis-equilibrium  with  and  without  the  presence  of  electro-magnetic  fields  have  been 
included  in  these  codes.  An  existing  parallel  code  environment  using  generalized  mesh  structure  (prisms, 
tetrahedra,  etc.)  developed  for  Computational  Electromagnetics,  was  adapted  to  model  MHD  with  (a) 
Finite  rate  chemistry  of  ionizing  air,  (b)  Thermal  Non-equilibrium  (Distribution  of  energy  across  various 
modes  of  storage),  and  (c)  Electromagnetic  effects.  The  effect  of  boundary  conditions  such  as 
conducting/non-conducting  walls,  applied  or  extracted  electric  power  and  so  forth,  render  this  predictive 
capability  invaluable  in  the  study  of  MHD  accelerators,  power  generators,  turbulence  control,  and 
integrated  analysis  of  fluid  mechanics  and  time-varying  electromagnetic  fields.  The  subject  of  inlet  flow 
control  for  hypersonic  vehicles  with  MHD  has  attracted  attention  due  to  the  possibility  of  flow 
modifications  without  the  need  for  moving  parts.  There  are  numerous  proposed  designs  to  reduce  drag, 
decrease  total  pressure  loss  and  enhance  air  mass  capture  using  MHD.  Additionally,  proposed  MHD 
based  concepts  can  even  serve  as  a  source  of  auxiliary  onboard  power,  and  as  means  to  increase  thrust 
produced  by  nozzles  in  hypersonic  vehicles.  The  code  development  activity  at  HyPerComp  has  resulted  in 
a  high  performance  toolkit  by  which  these  ideas  can  be  tested  computationally. 


Figure  1:  Unstructured  mesh  simulations  of  hypersonic  inlet 
with  MHD  interactions.  Above  left:  B  =  0,  below  left:  Bx  =  1 
Tesla,  cr  =  40  /Q/m  in  a  Mach  8  inlet.  Electrical  conductivity  is 
a  simple  temperature  dependant  function. 
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2.  Introduction 

The  application  of  plasma  physics  and  Mi  l D  to  aeronautics  has  recently  captured  public  attention  [1]. 
This  represents  a  resurgence  of  interest  from  the  initial  studies  in  hypersonic  MHD  [2]  and  the  pioneering 
work  in  the  then  Soviet  Union  ([3]  among  others,)  which  showed  shock  wave  attenuation  and  dispersion 
in  the  presence  of  plasma.  MHD  has  been  known  to  produce  numerous  beneficial  effects  in  hypersonic 
aerodynamics,  starting  from  enhanced  inlet  mass  capture,  drag  reduction,  better  combustion,  to  auxiliary 
power  generation  and  flow  acceleration.  Many  research  groups  have  reproduced  plasma-induced 
phenomena  such  as  shock  wave  attenuation  and  dispersion,  even  though  the  mechanisms  responsible  for 
these  effects  arc  not  clearly  understood  and  are  still  being  debated.  The  effect  of  cold  plasma  on  a  moving 
sphere  in  a  ballistic  range  is  shown  in  Fig.  [2].  The  Schlieren  photographs  here  show  gas  without  and  with 
a  glow  discharge.  The  attenuation  of  the  bow  shock  in  the  frame  with  the  plasma  turned  on  is  one  among 
numerous  illustrations  of  this  phenomenon.  Current  test  activities  within  the  United  States  (at  AEDC, 
John  Hopkins  University)  have  produced  drag  comparison  data  for  blunt  body  flows  with  and  without 
plasma.  These  studies  have  repeatedly  shown  an  increase  in  the  shock  stand-off  distance  and  a  decrease 
in  drag  in  the  presence  of  plasmas.  Tests  with  wind  tunnel  scaled  models  of  an  F-15  with  plasma  jet 
nozzles  in  the  nose  cone  have  been  conducted  at  IVTAN  (Russia).  A  series  of  flight  tests  (as  in  Fig.  [3]) 
have  been  proposed  by  Boeing  and  NASA-Dryden. 


Figure  3:  Illustration  of  F-15  equipped 
with  plasma  jet  nozzles  for  flight  test 

A  flurry  of  activity  has  also  emerged  in  the  area  of  hypersonic  propulsion  based  on 
magnetohydrodynamic  effects  which  has  returned  to  vogue  after  several  decades  [2,  4].  Herein,  a  Russian 
vehicle  concept,  referred  to  as  AJAX  (discussed  in  [5]),  will  bypass  a  portion  of  the  captured  airflow  and 
use  it  in  MHD  power  generation  and  subsequent  flow  acceleration.  The  MHD  generator-accelerator 
system  essentially  reduces  the  entropy  rise  of  air  diffusion  and  combustion  which  would  be  present  in  the 
unmodified  engine  flowpath.  Enhancements  to  this  vehicle  such  as  drag  reduction  by  injecting  plasma  at 
leading  edges  of  the  airframe,  and  increase  in  combustion  efficiency  by  plasma  interactions  have  been 
conceived  and  have  firmly  integrated  the  pursuit  of  the  following  two  issues  to  one  common  goal: 

(a)  Plasma  effects  in  subsonic  and  supersonic  flows  (vehicle  drag,  combustion  efficiency,  radar  cross 
sections,  etc.) 

(b)  The  use  of  MHD  in  high  enthalpy  flows  with  sufficient  electrical  conductivity  to  improve  inlet  mass 
capture,  shear  layer  mixing,  auxiliary  power  generation  and  nozzle  flow  acceleration. 

Thus,  there  is  a  strong  need  to  understand  the  various  physical  phenomena  and  their  impact  on  the  success 
of  the  hypersonic  propulsion  program.  However,  the  state  of  the  art  in  modeling  of  engineering  plasmas 
for  the  above  applications  is  not  adequate  to  comprehensively  include  all  the  relevant  physical 
phenomena.  Numerous  simplifying  models  exist,  which  addresses  some  of  the  physics  relevant  to 
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Figure  2:  Flow  past  a  sphere  without 
and  with  a  £low  discharge  plasma 
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hypersonic  MHD/propulsion.  Examples  of  such  simplified  approaches  include  the  simulation  of  finite  rate 
chemistry  for  multidimensional  problems  in  complex  geometries,  solution  of  the  integrated  Navier-Stokes 
and  Maxwell  equations  including  the  effects  of  finite  electrical  conductivity  (for  low  Hartmann  numbers 
at  least,)  and  the  validation  of  elaborate  experimental  data  describing  the  relaxation  processes  taking  place 
in  high  temperature  air.  More  and  more  complex  linkages  between  physical  phenomena  are  continually 
being  modeled:  for  example,  coupling  of  discharge  physics  with  surface  effects  of  partially  conducting 
objects  [5,  6],  e-beam  interaction  with  vibrational  energy  modes  of  air  [7],  influence  of  thermal  layers  on 
shock  propagation  [8]  and  so  forth.  It  would  appear  by  all  indications  that  the  situation  has  arisen  when 
the  next  generation  of  CFD  codes  based  on  sophisticated  physical  models  is  both  urgently  needed  and 
within  reach.  Such  a  CFD  code  would  be  able  to  model  supersonic,  viscous  flows  in  complex  geometries 
with  spatial  and  temporal  variation  of  transport  properties  and  include  effects  of  finite-rate  kinetics, 
thermal  non-equilibrium  and  electromagnetic  effects.  Such  a  tool  would  be  invaluable  in  the 
design/optimization  of  a  hypersonic  vehicle. 

The  objective  of  this  research  has  been  to  develop  a  comprehensive  design/optimization  tool  for 
hypersonic  propulsion.  A  computational  tool  which  includes  state-of-the-art  models  describing  the 
various  physical  phenomena  (chemical/plasma  kinetics,  thermal  non-equilibrium,  energy  transfer 
mechanisms  and  electromagnetic  effects)  was  developed  during  the  course  of  this  research.  The  code  was 
developed  to  study  3-D,  unsteady  flows  on  unstructured  grids.  A  general  purpose  “pre-processor”  was 
developed  to  generate  the  Jacobian  and  source  terms  for  an  arbitrary  set  of  chemical  reactions.  This 
capability  enables  the  study  of  various  chemical  kinetics  models  without  any  code  modification.  A 
Poisson  solver  generates  the  Lorentz  force  and  Joule  heating  terms,  which  are  used  in  the  momentum  and 
energy  equations.  Internal  mode  dis-equilibrium,  which  can  be  of  great  significance  in  gas-discharge 
phenomena,  can  also  be  modeled  with  this  code.  The  general  architecture  of  the  code  allows  it  to  be 
extended  easily  to  include  models  (transport  properties,  chemical  kinetics  etc.)  with  varying  degrees  of 
complexity.  This  capability  enhances  its  utility  as  a  research  as  well  as  a  design  tool. 

3.  Physical  Phenomena  of  interest 

Several  physical  phenomena  are  of  importance  in  hypersonic  MHD  flows.  Generation  and  sustenance  of 
charged  particles,  interaction  of  these  charged  particles  with  magnetic  fields,  non-equilibrium  energy 
transfer  mechanisms,  transport  properties,  viscous  and  turbulence  effects  and  wall/electrode/sheath 
phenomena  are  some  of  the  important  physical  effects  of  interest  in  such  flows.  The  success  of  a 
hypersonic  propulsion  concept  largely  depends  on  the  efficient  generation  and  sustenance  of  charged 
particles.  This  is  due  to  the  fact  that  Lorentz  forces  generated  by  the  applied  magnetic  field  depends  on 
the  electrical  conductivity  of  the  flowing  gas  stream.  All  the  above-mentioned  effects  such  as  chemical 
kinetics,  non-equilibrium  energy  transfer  mechanisms,  wall/electrode/sheath  effects,  transport  and  viscous 
effects  have  an  impact  on  the  spatial  and  temporal  variation  of  electrical  conductivity.  Hence  a  realistic 
assessment  of  a  hypersonic  propulsion  design  concept  would  require  a  thorough  understanding  of  the 
above-mentioned  phenomena  and  their  interaction.  Physical  phenomena  of  importance  in  the  study  of 
hypersonic  MHD  are  discussed  next. 

3.1  Mechanisms  to  create  and  sustain  ionization: 

A  practical  source  of  ionization  to  provide  a  large  enough  electrical  conductivity  is  essential  for  the  flow 
to  exhibit  MHD  interactions.  Hence,  a  key  issue  in  hypersonic  MHD  applications  is  ionizing  the 
incoming  air  stream.  The  requirements  for  producing  this  ionized  air  stream  are  very  stringent.  Large 
volumes  (~  1  m3),  low  gas  temperatures  (1000  -  2000  K),  high  electron  density  (~  1013  /cm3),  long 
residence  times  (>  10  msec),  low  power  budget  (1-10  MW/m3),  strong  magnetic  fields  (~  10  Tesla)  are 
the  main  requirements  for  these  air  plasmas.  Several  techniques  have  been  proposed  in  recent  literature, 
to  achieve  these  conditions.  In  the  past,  air  seeded  with  an  easily  ionizable  alkali  salt  [9]  has  been 
considered  alongside  the  use  of  electron  beams  [10].  Additionally,  high  voltage  nanosecond  pulses,  DC 
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&  RF  discharges,  E-beams  in  vibrationally  excited  gases  have  also  been  proposed  as  possible  means  of 
efficient  charge  particle  production.  Seeded  air  and  high  energy  E-beams  are  considered  to  be  viable 
means  of  charged  particle  generation  from  an  engineering  standpoint,  and  hence  were  studied  in  detail  in 
this  research. 

3.1.1  Seeded  air:  Seeded  air  flows  may  be  well  approximated  by  the  local  chemical  equilibrium 
assumption,  while  the  electron  species  may  still  be  in  thermal  nonequilibrium,  depending  upon  the 
incident  power  level  and  flow  properties.  Traditionally,  MHD  power  generation  concepts  have  involved 
seeding  combustion  effluents  with  some  easily  ionizable  salt,  such  as  Potassium  Carbonate,  Cesium, 
Sodium-Potassium  solutions,  and  so  forth.  There  are  many  practical  issues  involved  in  the  seeding 
process,  involving  multi-phase  flow,  deposition  of  seed  particles,  conglomeration  of  seed  particles  near 
cold  walls,  and  the  resulting  electrode  shorting,  to  name  a  few.  Nevertheless,  well-seeded  flows  exhibit 
high  electrical  conductivity  at  a  wide  range  of  pressures  and  temperatures.  The  seed  material  ionizes 
essentially  in  an  equilibrium  manner.  The  electrons  thus  produced  can  still  be  in  thermal  non-equilibrium 
with  the  rest  of  the  flow.  Predicted  values  of  electrical  conductivity  in  the  region  of  100  -  200  mho/m  can 
be  derived  at  temperatures  of  3000  K  and  flow  pressures  between  1  and  10  atm.  Munipalli  et  al.  [11] 
present  a  synthesis  of  a  seeded  air  calculation  procedure  which  uses  the  NASA-Lewis  CEA  code  to  obtain 
thermodynamic  properties.  Particle  laden  flow  is  simulated  including  the  effects  of  gravity,  and  seed 
particle  vaporization  (sublimation)  using  a  simple  linear  heat  transfer  analysis.  Highly  non-uniform  seed 
distributions  are  common,  greatly  reducing  effective  conductivity  of  gas  for  a  given  mass  of  injected  seed. 
Further,  losses  stem  from  relaxation  processes  resulting  in  heating  and  dissociation  of  seed  material.  Our 
present  thinking  is  that  seeded  air  flows  are  unusable  in  practical  aerodynamics,  except  for  ground  based 
testing  or  power  generation  applications  where  it  is  convenient  to  regain  some  of  the  injected  seed 
particles. 

3.1.2  e-beams:  The  modeling  of  incident  radiation  and  directed  energy  addition  has  been  proposed  at 
various  levels  of  complexity  in  the  recent  past.  This  could  also  include  the  use  of  pulsed  or  steady  RF 
(microwave)  discharge,  since  the  energy  transfer  mechanisms  still  involve  the  electron  species 
predominantly.  The  set  of  transport  equations  for  the  charged  particle  number  densities  as  proposed  by 
Macheret  [10]  to  evaluate  the  electrical  conductivity  a  for  such  situations  are: 

— -  +  V  •  Tc  =drJ  +  q i  +  kdNn_  -vane  -  §n+ne 
dt  11 

—■  +  V  •  f+  =  a|fe  1  +  4/  -  P„  ».  »+  -  P n+ne 
ot 

3h 

— ^  +  V  •  f_  =  -kdNn_  +Vane-  p/(  n_  n+ 

Ot 

where  the  fluxes  f  of  the  electrons,  positive  and  negative  ions  are  given  by  a  product  of  their  number 
densities  n,  and  their  effective  velocities.  The  gas  is  assumed  to  be  locally  neutral,  thus  implying 
n+=ne+n_.  (An  almost  identical  treatment  has  been  successfully  applied  to  model  an  RF  discharge  in 
[6]).  The  various  quantities  on  the  right  hand  side  of  this  equation  represent  the  ionization  due  to  plasma 
electrons  (a ),  e-beam  induced  ionization  rate  q„  collision^  detachment  of  electrons  from  negative  ions 
(k),  and  electron-ion  and  ion-ion  recombination  ( p  and  p„).  N  is  the  gas  number  density  and  v„  is  the 
electron-molecule  collision  frequency.  The  number  densities  may  then  be  used  in  any  model  describing 
the  electrical  conductivity  of  partially  ionized  gas.  The  simple  linear  model  proposed  in  [10]  relates  the 
number  density  computed  above  with  the  electron  charge  and  the  electron  and  ion  Hall  parameters  to 
obtain  the  electrical  conductivity  as: 
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It  is  desirable  to  use  an  expression  of  the  following  form  (derived  for  a  mixture  of  gases)  when  detailed 
species  chemistry  and  collisional  cross  section  data  is  available: 
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The  e-beam  induced  ionization  term  -  qi  -  is  the  critical  parameter  which  must  be  evaluated  to  obtain  the 
e-beam  power  balance  correctly.  This  quantity  is  computed  by  assuming  that  the  energy  deposition  profile 
(Qb )  is  close  to  being  Gaussian.  The  beam  induced  ionization  rate  is  then  given  by: 


where  the  denominator  represents  the  energy  cost  of  ionization  by  a  high  energy  beam  (in  this  case,  low 
energy  e-beams  are  known  to  be  highly  ineffective  and  lose  almost  all  their  energy  to  vibrational  modes 
in  the  first  few  milliseconds  upon  incidence.)  The  following  simple  expression  is  currently  being  used  in 
our  codes  to  model  the  electrical  conductivity  from  e-beams: 


1 

jbpYE„ ) 
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where  e  is  the  electron  charge,  m  is  the  electron  mass,  n  is  the  neutral  species  number  density,  kc  is  an 
electron  scattering  rate  constant,  kdr  (=  1.5e-7  cmA3/s)  is  a  dissociative  recombination  rate  constant,  Wi  is 
the  ionization  constant  (  =  34  eV/electron-ion  pair),  jb  (about  0.1  A/cmA2)  is  the  current  density  in  the 
electron  beam,  Eb  is  the  energy  of  the  electrons,  p  is  the  flow  density  and  Y  is  the  electron  stopping 
power.  The  energy  of  the  e-beams  can  be  varied  in  the  range  of  10  to  200  keV,  to  provide  high  enough 
electrical  conductivity  at  moderate  flow  static  pressures.  The  efficiency  of  electron  beams  diminishes 
with  increasing  pressure  due  to  the  inherent  nature  of  the  electrons  to  recombine  with  greater  ease  in  such 
situations. 

As  a  part  of  this  research,  an  E-beam  chemistry  model  proposed  in  [12]  was  implemented  to  study  the 
electrical  conductivity  and  MHD  effects  in  an  inlet  geometry  designed  for  Mach  8  operation  [13].  It  was 
concluded  from  this  study  that  an  E-beam  strength  of  1.5E24  m-3/s  (based  on  a  beam  current  density  of  5 
mA/cm2  and  beam  energy  of  30  keV),  yields  an  electrical  conductivity  of  about  5  mho/m.  This  low  value 
of  electrical  conductivity  is  inadequate  for  any  appreciable  MHD  effect.  The  beam  current  would  have  to 
be  increased  by  two  orders  of  magnitude  to  get  a  four- fold  increase  in  electrical  conductivity  [13],  hence 
making  it  infeasible  from  an  engineering  standpoint. 


3.2  Nonequilibrium  energy  transfer: 

3.2.1  Thermal  non-equilibrium:  The  preferential  redistribution  of  internal  energy  in  a  gas  to  vibrational, 
rotational,  electronic  and  translational  modes  may  be  modeled  with  the  knowledge  of  relaxation  laws,  and 
the  kinetics  of  dissociation  and  recombination.  Literature  is  abundant  with  data  pertaining  to  these 
phenomena.  Agreement  is  gradually  emerging  in  the  relative  importance  and  the  dominant  mechanisms 
which  are  responsible  for  global  phenomena  in  nonequilibrium  plasma  physics.  The  code  developed 
during  the  course  of  this  research  has  multi-temperature  capabilities,  namely,  the  translational-rotation 
temperature  (T),  vibrational  temperature  (Tv)  and  electron  temperature  (Te).  A  two-temperature  model 
[14]  provides  a  fair  representation  of  MHD  physics  in  core  flows  for  high  temperature  applications  (4000 
-  11000K),  at  which  range,  the  e-V  relaxation  time  is  on  the  order  of  nanoseconds  [15].  A  three- 
temperature  model  is  required  for  regions  with  lower  temperatures  or  in  regions  close  to  the  electrodes, 
where  the  e-V  relaxation  time  can  be  on  the  order  of  milliseconds  or  higher. 
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3.2.2:  Internal  mode  dis-eauilibrium  CState-Specific  Kinetics'):  Under  conditions  of  extreme  dis¬ 
equilibrium,  it  might  be  necessary  to  study  the  internal  modes  of  diatomic  gases  in  hypersonic  MHD. 
Under  such  situations,  the  internal  modes  (vibrational  states,  mostly)  need  to  be  modeled  as  individual 
species.  Processes  such  as  V-T  (Vibration-Translation),  V-V  (Vibration-Vibration),  Spontaneous 
Radiative  Decay  (SRD)  and  chemical  kinetics  must  be  included  in  the  source  term  of  these  master 
equations  to  model  the  production-depletion  of  the  internal  modes.  Recent  studies  have  shown  [16-18] 
that  associative  ionization  mechanism  in  an  optically  pumped  mixture  can  lead  to  significant  increase  in 
the  ionization  levels.  Vibrational  excitation  of  diatomic  air  species  have  a  profound  impact  on  the  electron 
removal  processes,  namely,  dissociative  recombination  and  attachment  of  02.  Vibrationally  induced 
detachment  of  electrons  from  02-  and  vibrationally  induced  heating  of  free  electrons  lead  to  enhanced 
ionization  levels  by  mitigating  the  effects  of  electron  attachment  and  electron-ion  recombination. 

The  MHD  code  architecture  developed  at  HyPerComp  has  the  framework  to  include  internal  mode  dis¬ 
equilibrium  (state-specific  kinetics)  coupled  with  flow,  chemistry  and  EM  fields.  The  vibrational  levels 
of  diatomic  species  would  be  modeled  as  individual  species  with  the  inclusion  of  V-T,  V-V,  SRD 
(spontaneous  radiative  decay)  terms.  Energy  transfer  into  the  higher  vibrational  manifolds  of  diatomic 
species  required  for  associative  ionization  and  mitigation  of  electron  attachment/electron-ion 
recombination  can  be  studied  with  this  model.  Such  a  study  would  elucidate  methods  to  control  the 
plasma  density  in  various  regions  of  the  flow.  Effects  of  diluent  species  on  mitigating  electron- 
attachment/recombination  can  also  be  investigated  using  this  detailed  state-specific  kinetic  model. 

3.2.3  Non-Maxwellian  Electron  Energies  -  Boltzmann  Solver  for  the  Electron  Energy  Distribution 
Function  fEEDFl:  Under  discharge  conditions,  low-pressure  molecular  plasmas  are  characterized  by 
several  complex  inelastic  processes.  Specifically,  the  energies  of  electrons  in  a  molecular  discharge 
couple  very  well  with  the  vibrational  states  of  tightly  bound  diatomic  species  via  inelastic  collisions. 
Consequently,  in  systems  where  V-V  exchange  collisions  dominate  V-T  exchanges,  the  vibrational  modes 
of  molecular  motion  are  out  of  equilibrium  with  the  translational  mode.  This  in  turn  causes  the  electron 
energies  to  also  be  out  of  equilibrium  with  the  translational  modes  of  the  gas  particles,  and  to  be  non- 
Boltzmann  without  a  clearly  definable  average  electron  energy.  Hence  it  is  incorrect  to  describe  the 
thermal  state  of  the  gas  using  an  electron  temperature  (which  assumes  a  Maxwellian  energy  distribution 
of  electron  energies,  described  by  the  electron  temperature).  The  electron  energy  distribution  and 
population  of  specific  excited  states  govern  the  rate  at  which  these  inelastic  collisional  processes  occur. 
Thus,  the  EEDF  determines  the  chemistry  and  hence  the  characteristics  of  the  plasma.  The  enhanced 
ionization  levels  reported  in  Ref.  [12]  is  perhaps  due  to  the  high-energy  tail  of  the  EEDF.  A  deeper 
understanding  of  flowing  plasmas  in  discharges  can  be  obtained  by  a  fully  coupled  approach  which 
involves  the  flow,  electromagnetics,  heavy  particle  energy  equation,  master  equations  governing  chemical 
kinetics  and  a  Boltzmann  equation  describing  the  EEDF. 


3.3  Electromagnetic  Fields: 

Electromagnetic  fields  play  an  important  role  in  the  study  of  hypersonic  propulsion.  The  Lorentz  body 
force  and  Ohmic  heating  terms  in  the  momentum  and  energy  equations  are  due  to  the  presence  of  the 
electromagnetic  fields.  Several  simplifying  models  have  been  used  to  compute  the  Lorentz  force  and 
Ohmic  heating  terms.  Analytical  and  computational  studies  on  MHD  in  the  literature  may  be  broadly 
divided  on  the  basis  of  the  magnetic  Reynolds  number  (Rem)  which  is  an  indication  of  the  relative 
magnitude  of  the  induced  magnetic  field  to  the  applied  field.  For  small  values  of  Rem  («  1)  the  flow 
conservation  laws  may  be  decoupled  from  the  Maxwell  relations.  This  is  sufficient  for  several  situations 
of  engineering  interest  (including  seeded  MHD  generators,  accelerators,  wind  tunnel  equipment  and  so 
forth.)  As  Rem  increases,  the  induced  field  becomes  progressively  stronger  and  can  have  a  significant 
impact  on  flow  physics.  Such  a  flow  permits  the  formation  of  complex  wave  systems  (Fast  and  slow 
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Magnteto-acoustic  waves  where  the  magnetic  field  and  flow  variables  may  jump  across  a  discontinuity, 
Alfven  waves  where  only  the  magnetic  field  rotates  in  a  direction  perpendicular  to  the  flow  and  other 
complex  interaction  patterns).  The  capture  of  these  wave  systems  places  stringent  requirements  on  the 
numerical  methods  employed.  Other  than  dissipative  and  dispersive  errors,  these  equations  suffer  from  a 
singular  eigenvalue  structure.  Powell  [19]  offered  an  empirical  fix  which  was  placed  on  firmer 
mathematical  foundation  by  Vinokur  [20]  which  has  been  used  successfully  in  upwind  schemes  in  recent 
times.  Some  of  the  approaches  of  coupling  electromagnetic  fields  with  the  flow  equations  are  discusses 
next. 

3.3.1  Source  term  approach:  The  most  straightforward  method  in  which  MHD  effects  may  be  introduced 
into  an  existing  CFD  code  is  by  the  addition  of  source  terms  corresponding  to  the  Lorentz  force  in  the 
momentum  equation,  and  Joule  heating  in  the  energy  equation  with  an  initially  prescribed  distribution  of 
electric  and  magnetic  field.  Such  an  approach  is  simple  to  implement  in  an  existing  fluid  flow  code.  The 
applicability  of  this  method  is  limited  by  the  size  of  the  magnetic  field  that  is  induced  in  the  flow  by 
virtue  of  the  currents  generated.  For  large  values  of  induced  fields,  the  magnetic  field  must  be  computed 
simultaneously  at  each  time  step.  Pressure  tends  to  be  underpredicted  by  such  a  source  term  based 
implementation. 


3.3.2  Constant  B.  with  spatially  varying  electric  field:  The  general  Ohm’s  law  including  an  electric  field, 
is  given  by  J  =  g(e  +  Vx  b).  For  situations  in  which  the  electric  field  may  be  deduced  from  a  potential 
(which  always  exists  when  the  magnetic  field  is  constant  in  time,)  may  be  rewritten  as: 

J  =  cr(-  Vip  +  Vx  B) 


In  order  to  determine  this  potential,  the  divergence  of  the  above  equation  which  is  zero,  by  the  law  of 
conservation  of  charge,  is  taken,  yielding  the  following  equation: 

V2<p  =V-(fxi?) 


Various  simplifications  may  now  be  applied  to  this  equation.  If  a  space  varying  magnetic  field  acting  in 
the  z-direction  is  applied  to  a  two-dimensional  velocity  field, 

V2  d(vB:)  S(uB.) 

V  dx  dv 


When  the  electrical  potential  is  computed,  it  is  readily  used  to  obtain  the  value  of  J  and  thus,  the  Lorentz 
force  JxB  may  be  deduced.  This  set  of  equations  must  be  augmented  by  the  appropriate  boundary 
conditions  described  in  the  following  section.  These  equations  may  be  written  out  trivially  for  a  general 
three-dimensional  situation  of  current  interest. 


3.3.3  Induced  magnetic  field  in  one  direction:  Several  authors  [21,22]  have  studied  situations  in  which  the 
applied  magnetic  field  is  oriented  along  a  fixed  direction.  This  leads  to  a  situation  wherein  the  induced 
electric  field  is  determined  by  Fleming’s  Left  Hand  Rule.  Consider  the  case  when  the  induced  field  Bx  is 
in  the  x-direction  and  the  applied  field  (0,By,Bz)  is  oriented  in  the  y-z  plane.  The  momentum  equation 
may  be  rewritten  in  terms  of  the  induced  field,  as  opposed  to  the  current  as  used  earlier.  This  and  the  time 
advancement  of  the  induced  field,  are  governed  by: 
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This  equation  without  the  time  derivatives  and  assuming  a  fully  developed  channel  flow  in  the  x-direction 
has  been  applied  with  insightful  results  in  [21,  22].  The  formulation  for  a  general  3-D  situation  where  all 
the  components  of  the  induced  field  exist  has  been  presented  in  Ref.  [23]. 
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3.3.4  Integrated  NS-Maxwell  equations:  A  full  complement  of  Navier-Stokes  and  Maxwell  equation  set 
has  been  solved  at  HyPerComp  for  problems  of  interest  in  compressible  flow  [24].  The  principal 
difficulty  in  solving  this  set  of  equations  is  to  maintain  a  divergence  free  magnetic  field  [25]. 
Characteristics  based  techniques  exist  in  compressible  flow  solvers,  designed  to  avoid  this  problem.  In 
incompressible  flows,  one  must  post-process  the  magnetic  field  to  remove  the  effects  of  non-zero 
divergence.  This  typically  involves  the  solution  to  a  Poisson-type  equation  and  may  be  computationally 
expensive.  However  the  equation  using  the  electric  potential  also  has  this  difficulty.  The  Maxwell 
equations  are  written  in  terms  of  the  magnetic  field  components. 

3.4  Near  wall  phenomena: 

The  effects  of  nonequilibrium  in  MHD  are  most  pronounced  near  the  electrode  region.  This  is  an 
important  region  from  the  standpoint  of  chemical  kinetics,  electromagnetics  and  thermal  non-equilibrium. 

3.4.1  Thermal  Non-Equilibrium  and  Chemical  Kinetics:  As  mentioned  above,  in  the  core  of  the  flow, 
vibrational  and  electron  modes  are  tightly  coupled  (for  temperatures  from  4000-1 1000K),  and  a  two 
temperature  model  proposed  by  Park  [14],  as  evidenced  from  the  work  of  Lee  [25]  may  be  adequate.  Near 
an  electrode  wall,  the  vibrationally  temperature  quickly  cools  to  the  boundary  value,  while  the  electron 
mode  loses  temperature  much  more  gradually,  well  into  the  sheath  region,  in  the  vicinity  of  a  few  mean 
free  paths  from  the  wall.  A  three-temperature  model  is  necessary  to  evaluate  the  electron  temperatures 
correctly  in  these  near-wall  regions.  Since  the  electron-induced  reaction  rates  depend  on  the  electron 
temperature,  an  accurate  prediction  of  charged-species  concentration  depends  on  the  correct  prediction  of 
electron  temperatures.  A  three-temperature  model  along  the  lines  of  [26]  has  been  developed  during  the 
course  of  this  work. 

3.4.2  Electromagnetic  Effects:  Sheaths  exist  in  the  near-electrode  regions.  Sheath  physics  in  hypersonic 
MHD  is  not  well  understood.  Sheath  thickness  is  on  the  order  of  a  few  Debye  lengths  and  can  be  varying 
from  a  few  microns  to  a  millimeter,  depending  on  the  electron  number  density  and  electron  temperature. 
Sheath  regions  are  characterized  by  high  electric  fields,  which  means,  highly  energetic  electrons  and 
hence  high  levels  of  ionization.  The  production  of  charged  species  in  the  sheath  region  has  a  profound 
impact  on  important  engineering  issues  such  as  electrode  erosion  and  electrode  shorting.  The  electrical 
conductivity  and  hence  Lorentz  forces  in  the  near-wall  region  impact  the  use  of  MHD  forces  in  flow 
control,  laminar  to  turbulence  transition  and  flow-separation. 

Wall  regions  also  present  a  challenge  in  terms  of  correctly  specifying  the  problem  from  a  numerical 
standpoint.  Boundary  conditions  have  to  be  correctly  specified  to  model  the  electromagnetic  effects 
accurately.  These  boundary  conditions  depend  on  the  choice  of  wall  material.  For  a  perfect  conductor, 
the  normal  component  of  the  magnetic  field  is  zero  and  the  jump  in  the  electric  field  has  no  tangential 
component.  For  a  perfect  insulator,  the  gradient  of  the  magnetic  field  is  zero  at  the  surface.  The  flow  of  a 
conducting  medium  in  the  presence  of  a  magnetic  field  produces  currents  that  must  form  closed  loops  or 
extend  indefinitely,  in  the  absence  of  charge  sources.  While  the  solid  walls  of  the  airframe  are  normally 
assumed  to  be  insulating  in  CFD  estimates  of  AJAX  performance,  even  mildly  conducting  airframe 
materials  cause  currents  to  leak  away  from  the  flow  (considering  the  very  low  achievable  values  of 
electrical  conductivity  of  air  in  flight,)  and  may  potentially  cause  major  changes  in  flow  compression  and 
shock  structures  on  this  account.  Hence  the  effect  of  conducting  walls  in  hypersonic  inlet  flow 
computations  must  be  considered  carefully.  Even  small  cracks  in  insulations  (4-5  orders  of  magnitude 
smaller  than  the  flow  length  scale)  have  large-scale  effects  on  the  flow,  and  the  influence  of  MHD  on 
flow  control  aspects.  We  have  initiated  studies  in  including  conducting  wall  phenomena  in  our  codes 
during  the  course  of  this  research  [27]. 
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3.5  Chemical  Kinetics:  A  chemical  kinetics  model  provides  critical  reactions  which  are  necessary  to 
describe  the  state  of  a  gas  under  a  given  set  of  thermodynamic  conditions.  Preference  is  usually  given  to 
simplified  reaction  mechanisms  which  involve  as  few  species  and  as  few  reactions  as  possible  which 
reproduce  the  correct  global  phenomena  such  as  ionization  fraction,  mixture  viscosity,  thermal 
conductivity,  vibrational  energy  absorption  rate  and  so  forth.  This  is  in  the  interest  of  saving 
computational  time  for  large  problems.  However,  these  simplified  chemistry  models  must  be  capable  of 
predicting  the  gas  species  composition  correctly.  These  models  must  be  capable  of  predicting  the  charged 
species  concentration  accurately,  in  order  to  get  realistic  estimates  of  the  electrical  conductivity.  An 
important  consideration  in  developing  chemical  kinetic  mechanisms  is  to  obtain  the  rate  constants  of 
individual  elementary  reactions.  Also  of  equal  importance  is  the  ability  to  obtain  thermodynamic  data 
(specific  heats  and  heats  of  formation)  of  the  individual  species  under  conditions  of  temperature  and 
pressure,  pertinent  to  the  operating  conditions  for  a  particular  problem.  Obtaining  reaction  rate  constants 
and  thermodynamic  data  for  high  temperature  low-pressure  conditions  (as  in  hypersonic  MHD)  is  a 
daunting  task.  Many  chemical  kinetic  models  have  emerged  over  the  years.  The  Park  model  (21 
reactions)  and  Dunn  and  Kang  model  (26  reactions)  [26]  have  been  found  to  provide  satisfactory  results 
for  a  wide  range  of  flight  conditions.  Using  ionization  mechanisms  such  as  e-beams  adds  further 
complications  and  demands  appropriate  kinetic  mechanisms  to  accommodate  them. 

There  are  several  numerical  challenges  in  chemical  kinetic  modeling.  Depending  on  the  reaction  rates,  a 
model  can  be  very  stiff  numerically,  thus  imposing  very  stringent  constraints  on  the  stability  and 
computational  time.  An  implicit  solver  provides  greater  stability  and  allows  the  use  of  large  time- 
constants.  A  major  impediment  in  the  use  of  an  implicit  scheme  involves  the  developing  the  Jacobian 
matrix  for  a  given  reaction  set.  This  process  can  be  considerably  time-consuming  and  also  very 
susceptible  to  errors.  To  alleviate  this  problem,  a  “pre-processor”  was  developed  to  automatically 
develop  the  Jacobian  and  source  terms  for  a  given  set  of  elementary  reactions. .  This  preprocessor  reads  in 
a  text  file  consisting  of  the  chemical  reactions  and  generates  the  Jacobian  and  source  terms.  The  Jacobian 
and  source  terms  can  be  included  in  the  solver  without  any  code  modification.  This  feature  greatly 
simplifies  the  process  of  testing  different  chemistry  models  and  is  thus  invaluable  in  obtaining  reduced 
chemical  kinetics  model  for  a  particular  set  of  conditions. 

3.6  Viscous  Flow  Features.  Turbulence  modeling 

Appreciable  electromagnetic  effects  will  occur  near  the  boundary  layer,  with  the  result  that  the 
temperature  profile  in  the  boundary  layer  will  change.  A  survey  of  turbulence  modeling  techniques  in 
MHD  may  be  found  in  Ref.  [28]. 

Turbulent  flows  comprise  of  vorticity  components  in  all  directions.  A  uniform  magnetic  field  acting  on  a 
turbulent  conducting  fluid  dampens  the  components  of  vorticity  perpendicular  to  itself.  A  scale 
independent  anisotropic  damping  of  the  velocity  field  is  caused,  which  makes  turbulent  eddies  essentially 
two-dimensional,  which  result  in  greatly  reduced  turbulence  and  are  damped  by  viscous  effects.  This  has 
been  confirmed  by  direct  numerical  simulation.  A  magnetic  field  acting  axially  in  a  pipe  with  turbulent 
flow  tends  to  damp  the  turbulence,  reduce  frictional  drag  for  a  given  flow  rate,  and  actually  raise  the 
critical  Reynolds  number  for  instability  of  laminar  flow. 

Enhancement  of  the  vorticity  in  the  flow  due  to  MHD  effect  and  modification  of  the  velocity  layer  near 
the  wall  due  to  interplay  between  change  in  the  temperature  profile  and  velocity  profile  can  be  the 
primary  source  of  change  in  the  turbulence  profile.  A  secondary  effect  may  be  the  change  due  to  electro¬ 
magnetic  jXB  term  in  the  momentum  equation,  which  will  change  the  turbulent  kinetic  energy  and 
turbulent  dissipation  equations  when  these  last  equations  are  obtained  by  taking  higher  moments  of  the 
momentum  equations.  While  turbulence  effects  might  be  important  in  hypersonic  MHD,  these  effects 
have  not  been  included  in  our  current  development. 
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4.  Governing  Equations  and  Numerical  Method 


4.1  Problem  Formulation: 

The  conservation  equations  governing  the  flow  of  weakly  ionized  gases  may  be  written  in  the  following 
form: 

—  +  V  •  (Fx,  Fy,  Fz)  +V  •  (Hx,  Hy,  Hz)  =  5 
dt 

in  which  the  fluxes  F  are  inviscid  or  convective,  while  H  are  viscous  or  dissipative.  The  source  term  S 
comprises  of  the  production  of,  and  interaction  between  the  various  flow  quantities.  This  equation  may  be 
expressed  in  the  integral  conservation  form  in  a  control  volume  O  bounded  by  80,  as  follows: 


+  (|  [F„ds  +  H„ds]=  ^SdO 
nan  n 

Fn  and  H„  in  the  above  equation  are  face  normal  flux  components.  The  vector  used  in  this  equation  is 
expanded  below.  An  index  s  is  used  for  the  various  gaseous  species  that  are  considered.  If  one  were  to 
consider  internal  mode  disequilibrium,  ps  would  comprise  of  various  internal  modes  (vibration/electronic 
levels)  of  individual  species.  We  show  the  three-temperature  formulation  along  the  lines  of  [26].  Energy 
addition  by  MHD  at  high  temperatures  can  bring  about  a  departure  from  equilibrium  between  the  various 
energy  modes.  There  are  three  energy  equations,  one  for  the  total  energy  conservation  and  the  other  two 
describing  the  vibrational  energy  of  the  diatomic  species  and  the  electron  energy.  The  following  matrices 
describe  the  conservation  laws: 
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The  transfer  of  energy  across  the  various  energy  levels  is  contained  in  the  non-equilibrium  energy 
production  rate  (6V  and  c6e .  The  source  term  in  the  vibrational  energy  equation  (cov)  contains  the  energy 

transfer  to  the  vibrational  mode  from  the  electronic  and  translational  modes  via  collisions,  and  vibrational 
energy  gained  or  lost  due  to  molecular  depletion  (dissociation)  or  production  (recombination).  The 
source  term  in  the  electron  energy  (c6g)  equation  contains  terms  due  to  the  electric  field  (Joule  heating) 

electron  pressure  gradient,  energy  exchange  with  the  translational  modes  of  heavy  particles  and  the 
vibrational  modes  due  to  collisions.  The  source  term  also  includes  energy  deposited  due  to  external 
radiations  and  energy  released  due  to  formation/depletion  of  electrons.  The  expressions  used  to  model 
these  features  may  be  found  in  [26].  The  three-temperature  formulation  is  explained  in  detail  in 
appendix-A.  In  this  work,  the  source  terms  in  the  electron  energy  equation  do  not  include  the  net  work 
done  due  to  the  electron  pressure  gradients.  This  is  justified  due  to  the  low  mass  fraction  of  electrons. 
Electron  velocity  is  assumed  to  be  equal  to  ion  velocity  and  hence  there  is  no  contribution  due  to  charge 
separation  effects.  Radiative  energy  addition  by  means  of  microwave  and  other  beams  can  possibly  be 
used  to  cause  an  increase  in  the  electron  energy  and  ionization.  This  is  modeled  by  the  terms  Qrad ■  Fluxes 
F„  and  H„  are  written  at  the  cell  face  in  each  control  volume.  Species  mass  fractions,  velocities  and 
temperature(s)  can  be  obtained  from  the  above  set  of  equations.  The  electromagnetic  effects  (current 
density,  electric  fields  and  potential)  are  obtained  by  the  solution  of  a  Poisson-type  equation  solved 
iteratively  at  each  time-step.  The  MHD  model  is  discussed  next. 

MHD  Model: 


In  this  work,  the  magnetic  Reynolds  number  is  assumed  to  be  small.  This  implies  that  the  induced 
magnetic  field  is  negligible  and  that  the  electric  field  may  be  derived  from  a  scalar  potential:  E  =  -V(p  . 
Ohm’s  law  may  then  be  used  to  relate  this  potential  to  the  current  density: 

/  =  a(- V(p  +  Fxs)-— (jxi) 

B 

Where  a  denotes  the  electrical  conductivity,  cox/B  is  the  Hall  parameter,  J  is  the  electric  current  density 
and  B  is  the  applied  magnetic  field  vector.  Since  the  induced  magnetic  field  is  neglected,  the  only 
electromagnetic  field  quantity  that  needs  to  be  computed  numerically  at  each  time  step  is  the  electric 
potential. 

When  the  divergence  of  the  above  expression  is  set  to  zero  and  is  re-written  in  terms  of  the  electric 
potential,  a  Poisson  type  equation  is  obtained,  which  must  be  solved  at  each  computational  time  step.  In 
the  results  presented  in  this  paper,  the  Hall-effect  terms  are  neglected  in  the  generalized  Ohm’s  law  given 
above.  Setting  the  divergence  of  the  current  density  (without  the  Hall  terms)  to  zero  leads  to  the 
following  expression: 

V  •  [cr  (-  Vq>  +V  x  i?)]  =  0 

Solution  of  the  above  Poisson  equation  yields  the  electric  potential  cp  and  electric  current  density  J.  The 
joule  heating  term  and  Lorentz  forces  in  the  energy  and  momentum  equations  can  thus  be  obtained. 


yPerComp,  Inc. 


L**d*r  in  High  Pdtfdttni int*  Computing 


AFOSR:  Inlet  performance 
enhancement  using  MHD 


Page  12  of  83 


4.2  Space/time  integration:  Numerical  Scheme 


A  point-implicit  finite  volume  scheme  is  used  to  solve  the  governing  equations  shown  earlier.  For  a  cell 
bearing  the  index  i,  volume  Q„  the  time  advancement  scheme  is  written  as: 


U‘+Al  =U‘  +CL~lr‘ 


Cl  - 

and 


where, 

i+—m'l 

n,  . 


+  Afvjsc  0.jMsrc 

The  right  hand  side  vector  r  is  the  flux  summation  based  on  quantities  at  time  level  t.  Inviscid  fluxes  are 
computed  from  a  Roe-type  procedure,  in  which  a  density  based  averaging  is  performed  similar  to  the 
perfect  gas  equivalent.  Viscous  and  dissipative  fluxes  are  computed  from  central  differences.  Details  of 
this  implementation  may  be  found  in  [15].  In  the  above  equation,  /  is  an  n  x  n  (n  is  the  total  number  of 
dependent  variables  in  the  vector  U)  identity  matrix  and  ML’  is  the  combined  Jacobian  including  inviscid, 
viscous  and  source  term  contributions. 


Reacting  and  plasma  flow  equations  are  characterized  by  widely  varying  spatial  and  time-scales.  The 
point-implicit  formulation  is  necessary  in  order  to  obtain  stable  solutions  for  a  stiff  system  of  governing 
equations.  The  formal  order  of  accuracy  of  this  numerical  procedure  can  be  enhanced  using  the  TVD 
approach,  as  described  by  Barth  [29],  using  local  gradients  computed  from  a  Gaussian  summation.  Work 
is  underway  to  develop  this  environment  into  a  higher  order  multi-physics  solver  using  the  Discontinuous 
Galerkin  technique. 


5.  Accomplishments 

5.1  Thermochemical  noneauilibrium 

Primarily  the  Park  two-temperature  model,  and  a  three-temperature  model  were  developed  in  the  course 
of  this  study.  In  the  two  temperature  model,  the  two  temperatures  used  were  T:  Translational  and 
Rotational  temperature,  Tv:  Vibrational,  electronic  and  electron  temperature.  In  the  three  temperature 
model,  the  three  temperatures  were  T:  (Translational  and  Rotational)  as  in  the  two-temperature  model, 
vibrational  temperature,  Tv  and  an  electron  temperature,  Te.  In  high  temperature  applications,  the  two- 
temperature  model  is  known  to  be  accurate.  The  energy  transmitted  to  the  electron  species  from  the 
electromagnetic  fields  is  quickly  transferred  to  the  vibrational  modes  and  at  a  much  slower  rate,  to  the 
translational  and  rotational  modes.  The  difficulty  arises  near  electrode  surfaces  or  situations  wherein  the 
temperatures  are  not  high  (as  in  E-beam  applications),  where  the  gas  temperature  might  be  less  than 
1000K.  An  example  of  the  applicability  of  the  two  temperature  model  is  a  high  temperature  MHD 
accelerator.  Figure  [5]  shows  the  axial  variation  of  temperature  for  a  two-temperature  model  in  an 
accelerator  channel.  It  is  seen  that  the  temperature  difference  between  the  vibrational  and  translational 
levels  remains  constant  through  the  length  of  the  channel.  This  difference  may  be  evaluated  from 
analytical  relations  as  in  Kerrebrock  [30],  thus  validating  the  nonequilibrium  model. 
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Figure  5  :  Nonequilibrium  Temperature  variation 
along  an  unseeded  air  MHD  accelerator  with  and 
without  the  MHD  fields 
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Figure  6:  Variation  of  electron  and  NO+  mole 
fractions  along  a  nonequilibrium  unseeded  air 
MHD  accelerator 


5.2  Equilibrium  models,  seeding 

The  NASA-Lewis  CEA  code  has  been  interfaced  with  a  Navier-Stokes  solver.  Ability  now  exists  to 
model  about  70  different  reacting  species,  with  about  1600  products,  including  multiple  phases.  As  a 
simple  illustration,  the  injection  of  liquid  droplets  of  Potassium  Carbonate  solution  in  water  has  been 
studied.  The  effect  of  finite  rate  evaporation  of  these  seed  particles,  and  their  settling  due  to  gravity  in  an 
MHD  channel  have  been  modeled.  The  results  from  this  study  were  originally  presented  in  Ref.  [1 1]. 

5.3  Integrated  NS-Maxwell  equation  modeling 

In  cases  where  the  Magnetic  Reynolds  number  is  large  (Rew  =<7|i0iC/oo  »  1,)  the  magnetic  field 
induced  in  the  flow  by  virtue  of  the  currents  generated  will  be  large,  and  a  system  of  combined  Navier- 
Stokes  and  Maxwell  equations  must  be  solved.  Recent  years  have  witnessed  a  significantly  enhanced 
understanding  of  the  computation  of  such  flows.  The  electrical  conductivity  values  encountered  in 
hypersonic  aerodynamics  are  rather  low,  and  a  study  was  performed  to  study  the  need  for  the  integrated 
NS-Maxwell  modeling  (ref.  [22]).  Three  different  numerical  methods  were  developed:  Roe’s  scheme, 
Xu’s  Kinetic  Flux  Vector  Splitting  scheme,  and  a  simple  central  difference  with  artificial  dissipation 
(Jameson-Schmidt-Turkel)  scheme.  With  varying  levels  of  difficulty  each  of  these  schemes  has  been  seen 
to  function  appropriately,  and  complex  MHD  wave  structures  have  been  resolved. 

5.4  Boundary  conditions 

Various  boundary  condition  procedures  have  been  developed  for  use  in  hypersonic  MHD. 

5.4.1  Effect  of  segmentation  of  electrodes:  In  applying  current  boundary  conditions  on  a  wall  comprising 
of  segmented  electrodes,  certain  simplifications  must  be  made.  If  the  computational  mesh  can  resolve 
each  electrode  sufficiently  (which  is  extremely  difficult  to  accomplish  in  real  problems,)  one  would  set 
the  normal  component  of  current  to  zero  in  the  insulating  region  and  the  tangential  gradient  of  the  electric 
potential  to  zero  on  a  conducting  surface.  However,  if  an  infinitely  segmented  electrode  arrangement  is 
assumed,  it  is  possible  to  deduce  the  following  condition  on  the  wall  [22]: 

=  _  n^P_ 
dy  dx 
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Induced  Magnetic  Field 


where  y  denotes  the  direction  normal  to  the  wall 
and  x  along  the  wall,  <p  is  the  electric  potential 
and  p  is  the  Hall  parameter.  This  expression 
seems  adequate  for  practical  applications. 

5.4.2  Magnetic  Field:  At  a  perfectly  conducting 
wall,  normal  component  of  the  magnetic 
induction  H  is  continuous  across  the  wall. 
Frequently  in  computations,  the  tangential 
component  of  this  field  is  set  to  zero,  which  is 
an  approximation  for  slender  bodies.  The 
magnetic  field  strength  inside  the  solid  body 
would  be  obtained  by  scaling  the  gas  magnetic 
field  strength  by  the  ratio  of  the  magnetic 
permeability  in  the  two  media.  The  induced 
magnetic  field  vectors  in  a  supersonic 
conducting  flow  past  a  ramp  is  shown  in  Fig. 
[7].  Here,  a  vertical  magnetic  field  is  applied, 
which  is  nearly  cancelled  tangential  to  the  wall  by 
tangential  induced  field  in  the  opposite  direction. 
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5.4.3  Various  treatments  of  induced  fields :  In  the  recent  past,  we  have  studied  in  considerable  depth,  the 
application  of  various  descriptions  of  MHD  flow  features.  Of  these,  the  main  ones  are:  (a)  The  induced 
electrical  field  approach,  (b)  Induced  B-field,  (c)  Induced  current  and  (d)  The  integrated  Navier-Stokes- 
Maxwell  equations.  A  few  sample  results  are  shown  here.  Fig.  [8]  shows  a  channel  flow  situation  in 
which  the  magnetic  field  is  ramped  up  from  0  to  0.35  Tesla  while  the  velocity  profile  flattens  and 
develops  an  M-shape  that  is  characteristic  of  such  flows.  Fig.  [9]  shows  the  effect  of  wall  thickness  and 
wall  conductivity  on  the  induced  magnetic  field  contours  in  a  fully  developed  channel  flow  in  3-D.  The 
development  of  the  3-D  M-shape  is  also  apparent  here. 


Figure  8  :  2-D  channel  flow  using  the  induced  electric  field  formulation.  Bz  is  zero  initially,  ramped  up 
to  0.35  Tesla  and  held  fixed  for  the  rest  of  the  channel.  U-profiles  are  shown  with  the  developing  M-shape. 
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Figure  9:  3-D  Channel  flow  using  the  induced  B-field  formulation.  Magnetic  field  has  a  bell  shaped 
peak  at  x  =  10.  “r”  is  the  ratio  of  wall  to  fluid  conductivities.  Plotted  alongside  is  the  velocity  profile 
along  the  cross  section. 


5.5  Wall  conduction  effects  in  Hypersonic  MHD: 

Existing  literature  in  the  area  of  hypersonic  flows  with  MHD  frequently  relies  on  the  assumption  that  the 
walls  bounding  ionized  flows  are  insulating,  or  have  some  pre-set  potential  distribution  that  determines 
the  current  within  the  fluid  medium.  When  the  walls  are  fully  conducting,  or  comprise  of  insulation  in 
which  imperfections  may  occur  (albeit  minute,)  literature  has  shown  numerous  situations  in  which  the 
computed  pressure  is  vastly  different  (by  more  than  an  order  of  magnitude,)  from  the  insulating  wall 
assumption.  In  airframes,  a  metallic  structure  is  more  of  a  norm  than  an  exception,  and  this  seems  to  call 
for  a  major  revision  of  MHD  performance  estimates  obtained  from  insulating  wall  approximations. 
Hence  a  part  of  this  research  was  focussed  on  understanding  the  impact  of  wall  physics  on  MHD  flows. 
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When  the  walls  are  thick,  computations  must  be  carried  out  within  the  wall  region  for  electromagnetic 
and  thermal  quantities  but  not  for  fluid  momentum  and  density.  This  involves  blocking  out  of  these 
regions  from  the  flow  solver.  An  index  array  denotes  the  material  type  in  each  cell.  Boundary  conditions 
are  applied  at  this  interface  rather  than  the  edges  of  the  computational  domain,  for  the  flow  solver.  When 
the  conducting  wall  surfaces  are  thick,  the  Boundary  Element  Method  (BEM)  may  be  used  to  model  the 
flow  of  currents  in  these  regions,  which  can  then  be  approximated  by  semi-infinite  boundaries.  In  the 
limit  of  thin  conducting  walls,  a  wall-tangential  Poisson-type  equation  may  be  coupled  with  the  flow 
solver,  thus  avoiding  the  need  to  generate  a  computational  mesh  inside  the  wall.  All  of  these  approaches 
are  being  presently  considered,  and  a  BEM  routine  is  being  developed  for  the  compressible  MHD  code 
environment.  Figure  [10]  shows  sample  computations  involving  conducting  walls.  Current  lines  in 
conducting  and  insulating  wall  situations  is  shown,  and  an  estimate  of  pressure  gradient  with  varying  wall 
conductance  ration  is  presented,  for  a  range  of  Hartmann  numbers  (these  estimates  are  based  on  the 
calculations  of  Tillack  [35]). 

5.6  Implementation  of  the  E-beam  chemistry  model: 

E-beams  are  considered  to  be  the  most  efficient  way  of  generating  charged  species  concentrations 
required  for  MHD  aerospace  applications.  A  major  thrust  in  this  research  was  to  develop  a  fully  coupled 
model  incorporating  the  effects  of  flow,  finite-rate  E-beam  induced  chemistry  and  electromagnetics.  The 
goal  was  to  compute  the  spatial  variation  of  electrical  conductivity,  based  on  charged-species 
concentration  and  use  it  with  the  MHD  model  to  obtain  the  Lorentz  forces  and  Ohmic  heating.  This  self- 
consistent  modeling  methodology  is  necessary  to  evaluate  the  potential  of  the  proposed  E-beam  concept 
for  several  MHD  based  applications  such  as  enhancement  of  inlet  mass  transfer,  drag  reduction  and  shock 
mitigation.  While  several  groups  have  studied  various  aspects  of  this  problem  [31-33],  none  of  these 
studies  have  attempted  to  do  so  in  a  fully-coupled,  self-consistent  manner.  During  the  course  of  this 
research,  self-consistent,  fully  coupled  model  was  used  to  study  a  3-D  Mach  8  inlet  geometry  under 
realistic  flight  conditions.  It  was  seen  that  for  values  of  beam  current  and  voltage  obtained  under  current 
engineering  constraints,  the  maximum  electrical  conductivity  is  too  low  to  achieve  MHD-based  flow 
control.  A  detailed  hypothetical  analysis  of  the  E-beam  model  (in  terms  of  beam  strength  and  extent  of 
the  E-beam)  was  conducted.  This  was  done  to  obtain  the  range  of  conditions  under  which  MHD-based 
control  would  be  possible.  Details  of  this  work  can  be  found  in  [13]. 

5.7  Implementation  of  the  three-temperature  model: 

A  simple  three-temperature  model  (details  given  in  Appendix)  was  also  developed  as  a  part  of  this  study. 
The  Jacobian  matrix  for  the  inviscid  flux  vector  F,  was  derived  along  the  lines  described  in  [26].  A 
simplified  2-D  geometiy  of  M-8  geometry  (as  in  Ref.  [13])  was  used  to  study  the  two-temperature  and 
three-temperature  models.  The  freestream  conditions  were  as  follows:  Tro  =  226. 51K,  PM  =  193.6  N/m2 
corresponding  to  conditions  at  30K  feet  (10K  meters).  Figure  1 1  show  the  variation  of  temperature  along 
the  surface  of  the  ramp  for  the  two-temperature  and  three-temperature  model.  It  is  seen  that  in  the  two- 
temperature  model,  the  vibrational  temperature  is  very  close  to  the  gas  temperature,  whereas  for  the  three- 
temperature  model,  the  electron  temperature  is  very  close  to  the  gas  temperature  while  the  vibrational 
temperature  is  virtually  unchanged.  At  low  temperature  and  pressures,  the  VT  relaxation  time  is  very 
large  (on  the  order  of  seconds)  according  to  the  Millikan- White  formula  [26],  however,  the  T-E  relaxation 
time  is  on  the  order  of  micro-seconds.  In  the  two-temperature  model,  since  the  electron  temperature  is 
assumed  to  be  equal  to  the  vibrational  temperature,  the  electrons  equilibrate  with  the  translational  mode. 
In  the  three-temperature  model,  however  it  is  clearly  seen  that  the  vibrational  temperature  is  low  (since 
the  energy  exchange  between  vibrational  mode  and  heavy  particles  is  on  the  order  of  seconds),  whereas 
the  electrons  equilibrate  with  the  translational  temperature.  This  comparison  shows  that  for  low 
temperature/pressure  cases  a  two-temperature  model  predicts  an  incorrect  vibrational  temperature. 
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T,  Tv  for  2-temp  model  T,  Tv,Te,  for  3-temp  model 

Figure  11:  Two  and  three  temperature  models  used  in  an  inlet  flow 

6.  Code  Architecture  and  Attributes 

The  development  of  the  HyPerComp  CFD/MHD  coupled  solver  is  based  on  the  code  architecture  of  an 
existing  unstructured  grid-based  CEM  solver.  This  finite-element  CEM  algorithm  has  been  implemented 
for  massively  parallel  and  scalar  computer  architectures.  Along  with  the  CEM  solver,  a  complete  suite  of 
user  interface  utilities  are  provided  in  the  computational  environment  to  go  from  a  CAD  model  to  the 
final  solution.  These  tools  are  shown  in  Figure  [12],  arranged  as  they  are  used  in  the  simulation  cycle. 

The  tools  include: 

•  The  Rockwell  Science  Center  developed  grid  generator  UNISG  for  performing  unstructured  grid 
generation. 

•  A  3-D  GUI  (graphical  user  interface)  for  pre  and  post-processing  called  UNSPREP.  This  tool  aids  the 
user  in  setting  up  grid-related  solver  parameters. 

•  A  help  sensitive  run  control  parameter  input  page  GUI. 

•  Domain  decompositioning  tools  for  parallel  simulations. 

•  The  solver,  which  has  been  ported  to  many  parallel  and  scalar  computer  architectures. 

•  A  collection  of  graphically  driven  post-processing  tools  for  performing  results  visualization  and 
solver  performance  analysis. 

The  software  written  by  HyPerComp  is  implemented  in  a  combination  of  C++,  C,  and  Fortran  and 
includes  C++  class  libraries  for  common  re-usable  functions.  It  should  be  noted  that  all  of  the  software  in 
the  solver  suite  is  either  freely  available  or  developed  in-house. 
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Tool  Sources: 

Blue  -  HyPerComp,  Inc.  developed  tools 
Black  -  Freely  available  external  tools 

Red  -  Commercial  tools 


Figure  12:  Integrated  electromagnetics  environment  at  HyPerComp,  Inc. 


7.  Summary  and  Recommendations  for  future  work 

7.1.  Summary  of  accomplishments: 

A  3-D  unstructured  code  capable  of  simulating  compressible,  supersonic  viscous  flows  has  been 
developed  during  the  course  of  this  research.  Effects  of  thermal  non-equilibrium  (multi-temperature 
models)  finite-rate  chemical  kinetics  and  electromagnetic  effects  have  been  coupled  with  the  flow  solver 
in  a  self-consistent  manner.  The  architecture  of  the  code,  derived  from  an  existing  CEM  solver,  makes  it 
amenable  to  be  ported  to  many  parallel  and  scalar  computers.  A  separate  preprocessor  is  used  to  generate 
the  Jacobian  and  source  terms  for  any  set  of  chemical  reactions.  This  preprocessor  reads  in  a  text  file 
consisting  of  the  chemical  reactions  and  generates  the  Jacobian  and  source  terms,  which  can  then  be 
included  in  the  solver.  This  code  framework  allows  the  user  to  study  any  arbitrary  set  of  chemical 
kinetics  (species  and  chemical  reactions)  without  the  need  for  any  code  modification.  Code  development 
for  studying  internal  mode  dis-equilibrium  (vibrational  non-equilibrium)  has  been  completed. 

This  research  has  resulted  in  the  first  study  of  a  validated  e-beam  chemistry  model  in  a  3-D  Mach  8  inlet 
geometry  under  realistic  operating  conditions.  Since  the  flow,  chemistry  and  electromagnetic  effects 
were  coupled  self-consistently,  realistic  estimates  of  electrical  conductivity  and  MHD  interactions  could 
be  obtained.  Several  important  engineering/design  issues  such  as  maximum  attainable  electrical 
conductivity  for  a  given  beam-strength  (beam  current  and  beam  energy  (keV)),  location  and  extent  of  the 
e-beam  were  investigated.  Based  on  engineering  constraints  on  achievable  beam  strengths  and  magnetic 
fields,  it  was  concluded  that  no  appreciable  MHD  interactions/flow  modification  could  be  achieved.  The 
charged  species  generation  by  E-beams  would  have  to  be  augmented  by  some  other  means  for  its 
successful  application  in  hypersonic  propulsion.  Vibrational  excitation  of  air  by  means  of  RF  or  laser 
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excitation  has  been  proposed  as  a  plausible  solution.  This  research  has  also  developed  the  first  3 
temperature  model  coupled  with  electromagnetic  effects  and  is  likely  to  lead  to  a  more  detailed  E-beam 
chemistry  model  with  rate  constants  based  on  electron  temperatures. 

7.2  Recommendations  for  future  research: 

Developments  to  the  code  can  be  made  in  the  following  areas,  namely  (a)  Physical  models,  (b)  Numerical 
techniques,. 

7.2.1  Physical  models:  While  the  current  code  addresses  several  physical  phenomena  relevant  to 
hypersonic  propulsion,  the  following  additions  is  likely  to  enable  a  better  understanding  of  the  underlying 
physics 

(a)  Turbulence  models:  Hypersonic  flowfields  in  the  fore-body  inlet  regions  of  proposed  NASP-type 
geometries  tend  to  be  dominated  by  turbulent  effects.  MHD  has  the  effect  of  suppressing  and 
laminarizing  turbulent  flows.  Two  equation  models  for  incompressible  MHD  turbulence  are  used 
widely  in  literature,  (see  e.g.,  [34]  for  a  description).  An  extension  of  these  models  to  compressible 
flows  has  been  undertaken  by  some  researchers  in  recent  times.  However,  this  problem,  like  much  of 
turbulence  modeling,  remains  open-ended,  particularly  due  to  the  lack  of  suitable  test  data.  In  an 
ongoing  research  effort  funded  by  the  Department  of  Energy,  we  are  investigating  closure  models  for 
incompressible  MHD  at  the  limit  of  large  Hartmann  numbers.  We  intend  to  extend  this  effort  into 
compressible  flows  in  due  course. 

(b)  Sheath  models:  Treatment  of  sheaths  in  hypersonic  MHD  can  be  a  daunting  task.  Since  our 
formulation  does  not  make  an  assumption  of  quasi-neutrality  in  the  computation  of  charged  species 
concentrations,  it  will  be  able  to  capture  sheath  effects.  However,  sheath-resolution  can  be  a  difficult 
problem.  Sheath  thickness  can  vary  from  a  few  microns  to  tens  of  microns  or  more,  depending  on 
several  factors  such  as  operating  conditions,  method  of  generating  charged  species  and  generation  and 
application  of  electromagnetic  fields.  A  series  of  parametric  studies  can  help  ascertain  the  sheath 
thickness  for  various  operating  conditions.  It  might  be  necessary  to  use  some  semi- 
empirical/analytica!  sheath  models  using  results  from  the  MHD  code  as  inputs. 

(c)  Internal  mode  dis-equilibrium  and  EEDF:  Code  development  for  the  inclusion  of  vibrational  dis¬ 
equilibrium  has  been  completed  during  the  course  of  this  research.  Knowledge  of  this  internal  mode 
dis-equilibrium  can  be  coupled  to  a  Boltzmann-solver  for  the  distribution  of  electron  energies.  The 
direct  coupling  between  the  electrons  and  internal  modes  can  be  studied  with  this  formulation.  The 
assumption  of  a  common  average  temperature  describing  the  electrons  can  thus  be  eliminated. 
Reaction  rate  constants  for  electron  induced  reactions  can  be  computed  on  the  basis  of  the  EEDF 
instead  of  the  Arrhenius  rates.  Important  physics  pertaining  to  reactions  induced  by  the  high-energy 
tail  of  the  EEDF  can  be  investigated  by  this  approach. 

(d)  Issues  related  to  relaxation  time  constants,  transport  properties  and  reaction  rate  constants:  There  are 
several  crucial  issues  related  to  accurately  predicting  the  electrical  conductivity,  relaxation  times  for 
VT,  V-E  and  T-E  processes,  transport  properties  such  as  thermal  conductivity  of  heavy  particles  and 
electrons,  viscosity  and  diffusion  co-efficients,  reaction  rate  constants  and  thermodynamic  properties. 
A  detailed  study  of  the  range  of  applicability  of  these  constants  has  to  be  ascertained.  This  is  of 
particular  significance  if  one  is  interested  in  pursuing  low-temperature  methods  of  producing 
ionization  such  as  E-beams,  RF  excitation  or  lasers. 

7.2.2:  Numerical  methods:  HyPerComp  is  currently  developing  higher  order  numerical  schemes  to  solve 
the  coupled  set  of  MHD  equations  using  the  induced  magnetic  field  formulation.  A  suite  of  codes, 
tentatively  christened  as  HOME  (Higher  Order  Multiphysics  Environment),  is  being  developed  based  on 
the  Discontinuous  Galerkin  approach  to  higher  order  numerics,  to  be  used  on  arbitrary  meshes,  with 
parallel  computers.  Higher  order  schemes  provide  the  ability  to  economically  resolve  complex  physics 
occurring  in  ordinarily  un-resolvable  length  scales,  such  as  Hartmann  layers,  acoustic-type  wave 
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propagation,  and  various  types  of  flow  instabilities.  It  is  of  imminent  interest  to  transfer  the  work 
developed  here  to  such  an  environment,  to  be  able  to  perform  full-airframe  computations  for  practical 
geometries. 
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8.  Interactions/Transitions 


Incompressible  MHD:  (Free  Surface  Liquids  for  Heat 
Removal)  In  this  DOE  funded  investigation,  HyPerComp  is 
studying  the  potential  for  using  electromagnetic  control  for 
heat  extraction  from  the  plasma  core  of  advanced  nuclear 
fusion  concepts.  Liquid  metals  (Lithium  or  Lithium  salts) 
flowing  along  the  walls  of  toroidal  magnetic  confinement 
fusion  devices  can  minimize  the  harmful  effects  of 
continued  exposure  to  neutron  bombardment  of  static 
walls,  and  at  the  same  time,  will  extract  heat  energy  more 
effectively  from  the  plasma.  An  added  benefit  is  the 
formation  of  Tritium,  which  is  used  regeneratively  as  fuel. 
A  critical  element  of  this  design  is  the  need  to  predict 
correctly  the  flow  of  liquid  metals  under  electromagnetic 
fields.  Three  dimensional  incompressible  flow  calculations 
involving  turbulent,  radiative  and  conductive  heat  transfer 
at  complex  free  surfaces  (involving  the  liquid  metal  and  an 
air/plasma  exterior),  will  be  the  ultimate  target  of  this 


Figure  1 1  :  Schematic  of  liquid  wall 
cooled  magnetically  contained  fusion 
device 


effort,  in  the  present  phase  and  beyond.  The  fusion  lab  at 

UCLA  is  the  immediate  beneficiary  of  this  research,  which  will  add  to  their  ongoing  DOE  funded  APEX 


(Advanced  Power  Extraction)  program.  A  code  environment,  tentatively  named  HIMAG  (HyPerComp 
Incompressible  MHD  solver  for  Arbitrary  Geometries),  has  been  developed  by  HyPerComp  under  these 
contracts.  HIMAG  is  already  making  significant  contributions  to  the  study  of  liquid  metal  flows  of 


interest  in  the  nuclear  fusion  community. 
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Appendix-A:  Three  temperature  formulation 

Details  of  the  three-temperature  model  used  in  our  work  are  described  below. 

Definitions: 

Translational  and  rotational  modes  of  heavy  particles  (T) 

Vibrational  energy  of  diatomic/polyatomic  molecules  (Tv) 

Translational  energy  of  electrons  (Tt.) 

Total  energy  equation/unit  mass 

ET  =  sum  of  translational  energies  of  heavy  particles,  kinetic  energy  due  to  flow,  vibrational  energy  of  all 
diatomics  and  electron  energy 

Ef  =  Y  ^"77^  +^Cses,y  +  CeCv,trEe 

s±e  z 

(Curve  fits  for  individual  species  account  for  the  heats  of  formation) 

Internal  energy/mass 
E  —  Y  cses  +  cecv  trTe 

s*e 

Vibrational  energy/mass 
Ev=YcseS'V 
Electron  energy/mass 

Ee  =cec%T 


The  total  pressure: 


n 


P  =  y!^RT  +p^rt 

M  M 

s*e  lvls  lvle 


Based  on  the  above  definition  for  pressure,  the  total  differential  for  pressure  can  be  written  as 


dp  =  RT±!!2*.  +  Rdr±%-  +  RT.%-  +  R%-', 


s*e  ^ S 


s#e  ^ $ 


Me  Me 


The  total  differential  for  pressure  can  be  written  in  terms  of  partial  differentials  as  follows: 

+^rdpE‘  +£dp‘ 

dp  E  dp  Ev  dp  Ee  dps 

Based  on  the  definition  of  internal,  vibration  and  electron  energies  given  above  the  differentials  for  the 
three  temperatures  can  be  written  as 
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=  dT 


ri 

dE  -  dEv  -  cedEe  -  £  (es  -  ev,s  )dcs 


" vJr 


ri 

dEv -^(eViS)dcs 


=  dT„ 


V,V 

&  =  dT ; 

e  e 


r 

^v,tr 

The  differentials  for  energies  in  terms  of  the  primitive  variables  can  then  be  written  as 

Jr,  (dpE  -  Edp  -  ( udpu  +  vdpv  +  wdpw)  +  ( u 2  +  v2  +  w2)dp) 

dE  =  — — - - 

p 

dE  _  (dpEv-Evdp) 

P 

dE  _  (dpEe-Eedp) 

P 

dc,  =  (dps~csdp  ).  =»  p dcs  +  csdp  =  dp s 


Substituting  the  above  expressions  for  dT,  dTe,  dE,  dEv,  dEe  in  the  total  differential  for  pressure,  we  have 
the  following  definitions. 

__dp_^  R  ^  P r 


P 

(P  = 

P 

y,= 


dpE  pCv  lr  s±e  Mr 

dp  =  RPe  qc 
dpEe  pCl,rMe  e 

dpEv 

dP  _^Tq  ,  p  (»2+V2+^2) 


dps  Ms 


P«,  +  P 


Where,  Tq  =  Te  if  5  =  e  and  Tq  -  T  if  s*e 

Based  on  these  definitions  the  Jacobian  of  the  inviscid  matrix  F  can  be  written  as: 
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Source  term  for  vibrational  energy  in  a  two-temperature  model: 


wv  =  i 


Y,wsDs~  Z 


s-mol 


r-elct.imp 


Rf,r-Rb,rVr 


Z  P-s 

s~mol 


E  v,s  Ev,s 


i's) 


+  3t>eR(T-Tv)Y, 


es 


s*e  M s 


+  farad 


Source  terms  for  the  vibrational  and  electron  energies  in  the  three-temperature  model  is  given 
below. 


We=3PeR{T-Te)Y^-  I  P5 


f  ** 

E  v,s  “  Ev,s 


s^e  M $  s-mol 
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The  definition  of  various  terms  in  the  above  equations  is  defined  below. 

^wsDs  =>  Average  vibrational  energy  per  unit  mass  created  or  destroyed 

J?,SRfs~Rb’-  X  Average  energy  per  unit  mass  created  or  destroyed  due  to 

change  in  electron  number  density 
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3p eR(r  -  rv)Z—  =>  Energy  exchange  between  heavy  particles  and  electrons 
s*eMs 

Qrad  -  A-V.w  =>  source  terms  due  to  external  radiation  and  electron  pressure  gradient 
f  **  \ 

K  '  =>  energy  exchange  between  electrons  and  vibrational  energy 
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Abstract.  We  present  here  our  recent  studies  in  the  computation  of  electrically  con¬ 
ducting  air  flow  in  hypersonic  inlets  in  the  presence  of  electromagnetic  fields.  In  nu¬ 
merical  simulations,  one  has  the  opportunity  to  study  such  flows  under  various  levels 
of  abstraction.  Models  studied  range  from  linear  Lorentz  force  terms  to  induced  field 
formulations  using  reduced  as  well  as  a  full  treatment  of  Maxwell  equations.  Real  gas 
phenomena  such  as  chemical  and  thermal  nonequilibrium  are  considered.  Electrical 
conductivity  enhancement  via  seeding  and  the  use  of  electron-beams  for  significant 
interactions  to  occur.  A  comparative  study  of  numerical  models  and  physical  issues  is 
presented. 


1  Introduction 

Interest  in  MHD  applications  to  aeronautics  dates  back  to  the  1950s  (e.g.,  Resler 
and  Sears  [9],  following  its  active  study  in  problems  related  to  interstellar  gas- 
dynamics,  and  ground  based  power  generation.  The  resurgence  of  interest  in  the 
present  times  is  owed  largely  to  the  Russian  AJAX  concept  in  which  a  scramjet 
inlet  flow  field  is  bypassed  through  an  auxiliary  MHD  power  generator,  reduc¬ 
ing  the  net  flow  compression  required.  Benefits  arise  from  lower  drag,  increased 
combustion  efficiency  via  plasma  injection,  generation  of  additional  power,  flow 
control  without  the  need  for  variable  geometry  in  the  inlet  as  well  as  the  nozzle. 
Ref.  [1]  may  be  consulted  for  a  recent  account  of  these  concepts. 

Significant  advances  have  been  made  in  the  development  of  numerical  schemes 
for  MHD  in  the  recent  past.  Powell  [8]  developed  a  method  which  provides  for 
the  proper  advancement  of  solenoidal  magnetic  fields  in  time,  thereby  enabling 
the  development  of  effective  Riemann  solvers  in  MHD.  DeStreck  et  al.  [2],  Linde 
[4],  among  others,  have  performed  high  accuracy  computations  and  attempted  to 
interface  current  numerical  methods  with  classical  characteristics  based  analysis 
of  MHD.  Direct  applications  to  aerospace  engineering  problems  have  been  taken 
up  by  several  researchers  (e.g.,  [3],  [5]). 

A  hierarchy  of  gas  models  and  MHD  models  have  been  considered  in  the 
present  work.  The  intent  is  to  develop  a  three  dimensional,  unstructured  (hy¬ 
brid)  mesh  code  to  run  on  parallel  machines  for  general  time  accurate  MHD 
problems.  An  existing  code  environment  developed  for  use  in  computational 
electromagnetics  is  being  modified  for  MHD.  Refs.  [6]  and  [7]  may  be  consulted 
for  further  details  on  past  activities.  Recent  results  pertaining  to  inlet  studies 
are  presented  here  with  a  brief  summary  of  physical  models  used. 
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2  Gas  Models 

An  important  distinguishing  feature  of  magneto-aerodynamics  as  encountered  in 
hypersonics  is  the  low  electrical  conductivity  and  the  perennial  quest  to  increase 
its  value.  The  choice  of  a  gas  model  is  closely  related  to  the  manner  in  which 
this  electrical  conductivity  is  enhanced  or  produced.  Among  the  most  popular 
options  to  increase  the  electrical  conductivity  are  :  (a)  Seeding  the  flow  with 
an  easily  ionizable  salt,  (b)  Using  electron  beams  or  other  forms  of  radiation  to 
increase  the  concentration  of  mobile  electrons  in  the  flow. 

(a)  Perfect  gas:  The  perfect  gas  assumption,  while  being  physically  inaccu- 
rate  in  general  MHD  flows  of  engineering  interest,  is  preferred  due  to  its  simplic¬ 
ity  and  speed  in  numerical  simulations.  Parameters  such  as  the  ratio  of  specific 
heats  (7)  and  the  gas  constant  can  be  varied,  to  suit  the  mean  flow  proper¬ 
ties.  Sutherland’s  law  type  of  temperature  dependence  of  viscosity  and  thermal 
conductivity  is  used. 

(b)  Chemical  equilibrium:  We  have,  in  the  past,  (Ref.  [6])  used  the  equilib¬ 
rium  air  assumption  in  order  to  study  high  temperature  air  seeded  with  Potas¬ 
sium  Carbonate  particles.  Seeding  provides  an  effective  mechanism  to  enhance 
flow  conductivity  at  temperatures  in  the  vicinity  of  3000  K.  Simple  corrections 
for  multiphase  effects  of  heat  transfer,  phase  change  and  particle  settling  due  to 
gravity  have  been  made  in  our  equilibrium  air  simulations.  A  set  of  routines  de¬ 
veloped  at  NASA-Glen  Research  Center  are  being  used  to  model  the  equilibrium 
composition  and  properties  of  high  temperature  air  with  and  without  seeding, 
when  finite  rate  effects  are  not  significant. 

(c)  Thermochemical  nonequilibrium:  Seed  particles  result  in  a  weight  penalty, 
contaminate  the  flow,  and  are  rarely  recoverable  for  cyclical  usage.  Currently, 
there  is  interest  in  ionizing  air  by  means  of  electron-beams  and  other  forms  of 
radiation.  Such  methods  of  ionization  tend  to  be  dominated  by  nonequilibrium 
processes.  Electron  collision  and  recombination  times  are  comparable  to  the  flow 
time  scales. 

Hypersonic  flows  at  high  altitudes  possess  their  own  charge  of  nonequilibrium 
features.  Energy  distribution  among  the  various  modes,  as  well  as  chemical  reac¬ 
tions  among  gaseous  species  take  place  at  a  finite  rate.  In  our  work,  we  have  used 
the  two-temperature  model  developed  by  Chul  Park.  In  this  mode,  the  transla¬ 
tional  and  rotational  temepratures  are  characterized  by  one  temperature  scale 
(say  T,)  and  the  vibrational  and  electronic  modes  are  characterized  by  Tv.  This 
model  degenerates  in  accuracy  near  electrodes,  where  the  electron  temperature 
differs  significantly  from  vibrational  and  translational  temperatures. 

Electrical  conductivity  from  e-beams  is  frequently  expressed  as: 


where  e  is  the  electronic  charge,  m  is  the  electron  mass,  n  is  the  neutral  species 
number  density,  kc  is  an  electron  scattering  rate  constant,  kar  (=1.5  e- 7  cm3/ s) 
is  a  dissociative  recombination  rate  constant,  is  the  ionization  constant  (=34 
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eV),  jb  (about  0.1  A/cm2)  is  the  current  density  in  the  electron  beam,  Eb  is 
the  energy  of  the  electrons,  p  is  the  flwo  density  and  Y  is  the  electron  stopping 
power.  The  energy  of  the  e-beams  in  the  range  of  10-200  keV  is  needed  to  provide 
high  enough  electrical  conductivity  at  moderate  flow  static  pressures. 

3  MHD  Models 

When  the  magnetic  Reynolds  number  (Rem  —  VooLapo)  in  a  flow  is  much  greater 
than  1,  sizeable  magnetic  fields  are  induced  in  the  flow  from  the  currents  that 
are  generated  by  virtue  of  its  movement  in  a  magnetic  field.  These  fields  are  not 
significant  in  aerodynamics,  where  the  electrical  conductivity  is  often  very  low. 

(a)  Inductionless  MHD:  In  the  limit  of  low  magnetic  Reynolds  number,  the 
Navier  Stokes  and  Maxwell  equations  are  decoupled.  A  shock  tube  problem  has 
been  simulated  with  MHD  included  a s  source  terms  where  a  magnetic  field  is  ap¬ 
plied  in  the  z-direction  and  no  electric  or  magnetic  field  is  induced.  As  remarked 
earlier,  the  magnetic  field  is  seen  to  retard  the  shock  waves  and  dissipate  the 
discontinuities.  A  constant  value  of  electrical  conductivity  (40  /Q/m)  has  been 
assumed.  While  the  codes  we  have  developed  include  Hall  effect,  consider  for 
the  present  Ohm’s  law  neglecting  Hall  effect  and  electron  pressure  gradient: 
j  =  <r(-V0+  VxB) 

This  current  must  be  divergence  free,  as  per  Maxwell’s  equations.  In  order  to 
satisfy  this,  it  is  customary  to  rewrite  the  above  equation  in  terms  of  an  electric 
potential  ( E  =  -V0)  and  take  the  divergence  of  both  sides.  The  following 
potential  relation  may  then  be  obtained  for  the  electric  potential.  Flows  with 
and  without  this  induced  potential  have  been  studied. 

V  •  j  =  0  =►  V2<£  =  V (log a)  *  (-' Vcj>  +  V  x  B)  4-  V  •  (V  x  B)  (1) 


Fig.  1.  Density  profile  in  a  shock  tube  flow  with  MHD  as  source  terms  (left)  and 
Pressure  profile  showing  the  effect  of  finite  a 
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(b)  The  induced  field  formulation:  In  the  general  case  of  large  electrical  con- 
ductivity,  induced  magnetic  fields  are  significant  and  present  a  wave  structure 
that  cannot  be  captured  by  the  previous  formulation.  Fig.  1  shows  a  shock  tube 
problem  with  a  set  to  zero,  a  finite  value  (1000  /O/m)  and  oo.  The  presence  of 
additional  waves  due  to  MHD  and  their  eventual  damping  due  to  finite  a  is  seen 
here. 


4  Numerical  Methods 


In  the  case  of  the  inductionless  approximation,  MHD  terms  are  included  as  source 
terms  in  the  Navier-Stokes  (NS)  equations,  and  as  constant  terms  in  the  Jaco¬ 
bian  matrix  when  an  implicit  scheme  is  used.  The  Poisson  equation  for  electric 
potential  is  solved  by  relaxation.  Heavy  under-relaxation  (u  =  0.01)  is  some¬ 
times  necessary  for  strong  shocks.  Roe’s  scheme  with  MUSCL  type  evaluation 
of  fluxes  is  used  to  obtain  second  order  accuracy  in  spatial  coordinates.  Where 
time  accuracy  is  desired,  a  second  order  Runge-Kutta  time  stepping  was  used. 

To  solve  the  coupled  NS-Maxwell  equations,  the  Powell  [8]  formulation  is 
used,  wherein  a  source  term  proportional  to  the  divergence  of  B  is  introduced 
on  the  right  hand  side.  Three  different  solution  procedures  were  used:  (1)  Xu’s 
Kinetic  Flux  Vector  Splitting  [10],  (2)  Roe’s  scheme  as  presented  by  Linde  [4], 
and  (3)  a  central  difference  scheme  with  artificial  dissipation  based  on  the  mag¬ 
netic  pressure.  Neumann  stability  of  the  finitely  conducting  MHD  equations 
depends  upon  two  parameters  v  and  r: 


v 


(u  +  cf) 


At 

Ax' 


r  = 


where  u  +  c/  is  the  fastest  wave  speed  at  a  cell  interface.  For  stability  of  an 
explicit  scheme: 

i/2  <  2r,  r  <  ~ 


5  Results  and  Discussion 

(a)  Mach  10  inlet  at  off-design  conditions:  A  two-ramp  inlet  designed  for  opera- 
tion  at  Mach  10  is  shown  in  Fig.  2  at  M  =T  10  and  6.  At  a  Mach  number  of  6,  not 
only  are  the  shock  strengths  too  large  for  the  compression  obtained,  but  also  the 
intake  mass  capture  is  greatly  reduced.  An  equilibrium  seeded  flow  study  and 
a  simulation  using  e-beams  to  enhance  electrical  conductivity  of  the  flow  are 
shown  in  Fig.  3.  When  a  magnetic  field  of  5  Tesla  is  applied  in  the  x-direction, 
the  mass  capture  is  enhanced  by  about  9%,  and  the  total  pressure  loss  is  reduced 
by  about  5%. 

(b)  Coupled  NS-Maxwell  simulations:  Ideal  MHD  simulations  of  flow  past  a 
5°  ramp  at  Mach  6  are  shown  in  Fig.  4.  A  magnetic  field  (of  different  strengths) 
is  applied  in  the  y-direction.  An  insulating  solid  wall  condition  is  imposed  on 
the  ramp  wall,  and  a  ratio  of  relative  magnetic  permeabilities  between  the  flow 
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Fig.  2.  2D  inlet  at  design  condition  M=10  (left)  and  off-design  at  M=6  (right) 


Fig.  3.  Seeded  flow  with  applied  Bx  (left)  and  Flow  ionized  by  e-beams  (right) 


and  the  solid  wall  was  taken  as  14.  The  continuity  of  magnetic  induction  vector 
H  —  fioB  was  imposed.  The  formation  of  fast  and  slow  magnetosonic  shock 
waves  is  seen  in  this  illustration.  The  two-dimensional  nature  of  this  simulation 
does  not  permit  the  formation  of  an  Alfven  wave.  At  low  values  of  a,  these 
discontinuities  merge  into  seemingly  a  continuous  (nearly  linear)  transition  of 
pressure  from  free  stream  to  the  wall  condition  after  the  leading  shock  wave. 

(c)  Nonequilibrium  simulations:  A  model  of  air  involving  11  species:  N>  0, 
N2 ,  02,  no,  iV+,  0+  JV2",  of,  NO+  and  and  two  temperatures  as  described 
earlier  is  used  in  our  CFD  studies.  21  chemical  reactions  between  these  species 
including  the  ionization  of  NO  (which  is  the  significant  source  of  electrons  in 
hypersonic  flows  considered,)  are  studied.  Refs.  [6]  and  [7]  present  more  details 
about  the  nonequilibrium  studies.  Here,  we  present  a  double  ramp  inlet  designed 
for  Mach  8,  operating  at  low  supersonic  (M  =  2.68)  speed.  A  magnetic  field 
projecting  out  of  the  plane  of  the  paper  is  applied  centered  at  x  =  4m,  spanning 
lm.  The  upper  and  lower  boundaries  are  assumed  to  be  non-conducting  walls. 
Current  direction  here  will  be  tangential  to  the  walls.  Electrical  conductivity 
was  enhanced  by  introducing  a  source  term  in  the  electron  continuity  equation 
that  is  active  in  the  region  where  the  magnetic  field  is  non-zero. 

Fig.  5  shows  the  variation  of  temperature  along  the  wall  for  cases  with  and 
without  MHD.  The  vibrational  temperature  absorbs  the  Joule  heating  first, 
and  rapidly  equilibriates  with  the  translational  temperature.  In  this  situation, 
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Fig.  4.  Effect  of  magnetic  field  on  magnetosonic  shocks 


nonequilibrium  processes  are  not  critical.  Electric  field  is  computed  from  the  po¬ 
tential  equation  as  discussed  earlier.  An  interesting  feature  in  this  flow,  is  that 
if  the  electric  field  is  not  computed,  the  magnetic  field  makes  the  flow  unsta¬ 
ble,  as  in  an  inlet  unstart.  A  strong  shock  is  formed,  down  stream  of  which  the 
flow  is  subsonic,  and  propagates  upstream  and  “unstarts”  the  inlet.  This  is  seen 
in  figs.  6.  The  solution  obtained  with  the  electric  field  is  stable  and  converges 
smoothly.  When  this  field  is  neglected,  it  appears  that  no  steady  state  solution 
is  possible. 

It  is  a  tedious  task  to  validate  MHD  flow  features,  due  to  their  highly  non- 
intuitive  nature.  A  systematic  study  is  underway  to  unite  various  research  teams 
to  set  up  test  problems  and  analyze  the  effects  of  physics  and  numerics  in  nu¬ 
merical  MHD.  A  study  of  this  nature  in  ideal  MHD  seems  to  be  already  under 
progress,  as  described  in  Ref. [2].  Hypersonic  MHD  is  dominated  by  additional 
physical  issues  whose  synergy  seems  to  be  yet  understood  and  exploited  for 
engineering  use. 
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Abstract 

In  this  paper  we  describe  the  development  of 
simulation  capabilities  to  perform  computations  of 
magneto-hydrodynamic  (MHD)  flow  in  large 
complex  domains.  Numerical  simulation  tools 
available  in  the  present  time  are  heavily  inadequate  to 
realistically  simulate  the  motion  of  conducting  gases 
in  the  presence  of  electromagnetic  fields.  Realism, 
here  refers  to  the  inclusion  of  finite  conductivity, 
thermal  nonequilibrium,  charge  seperation,  chemical 
reactions,  dissipative  effects  and  the  ability  to 
simulate  large  induced  magnetic  fields.  Much 
understanding  has  been  gained  in  the  numerical  study 
of  ideal  fluid  MHD  in  recent  times.  Likewise,  real 
gas  computations  involving  thermochemical 
nonequilibrium  and  gas-kinetic  modeling  of  transport 
properties  are  also  plentiful  in  recent  literature.  This 
work  attempts  to  combine  these  disciplines  to 
develop  a  comprehensive  numerical  treatment  of 
MHD.  This  paper  is  organized  into  two  parts.  The 
first  deals  with  an  assessment  of  physical  models  and 
available  data  in  MHD  modeling.  The  second 
presents  CFD  techniques  used  to  model  the  integrated 
Navier-Stokes  -  Maxwell  equation  set  in  multiple 
dimensions. 

Physical  Phenomena  in  Real  Gas  MHD 

Generation  and  sustenance  of  a 
The  operation  of  any  device  based  on  MHD  depends 
critically  upon  the  level  of  electrical  conductivity  that 
may  be  attained.  Electrical  conductivity  may  be 
enhanced  principally  by  the  following  mechanisms: 

1.  Increase  in  gas  temperature:  The  bulk  of  the 
electric  current  in  a  conductor  is  carried  by  the 
electrons  due  to  their  low  mass.  Raising  the 
temperature  of  a  gas  causes  an  increase  in  the 
degree  of  ionization,  thereby  increasing  the 
electron  number  density.  This  form  of  ionization 
is  frequently  referred  to  as  equilibrium 
ionization.  The  electron  concentration  results 
predominantly  from  the  ionization  of  Nitric 
Oxide  (NO)  molecules,  forming  NO  ions,  until 
about  6000  K.  The  usable  values  in  the  vicinity 
of  100  mho/m  are  attained  at  prohibitively  high 


values  of  gas  temperature  (over  7,000  K).  The 
need  is  therefore  felt  for  additional  mechanisms 
to  enhance  electrical  conductivity. 

2.  Nonequilibrium  ionization  by  preferential 
elevation  of  electron  energy:  Diatomic  Nitrogen 
remains  the  dominant  diatomic  species  carrying 
the  vibrational  energy  over  much  of  this 
temperature  range.  Energy  exchange 
mechanisms  between  Nitrogen  and  the  electron 
species  have  been  the  subject  of  several 
investigations  spanning  over  half  a  century.  It  is 
fairly  certain  at  the  present  time  that  under 
moderate  electromagnetic  fields,  the  electron 
energy  is  in  close  equilibrium  with  the 
vibrational  energy  of  air,  when  compared  with 
other  forms  of  energy  storage,  e.g.,  rotational  and 
heavy  particle  translational  modes.  In  air,  energy 
can  be  added  preferentially  to  the  electron 
species  by  means  of  Joule  heating.  An  electric 
current  with  a  component  parallel  to  an  applied 
electric  field  heats  the  gas  at  a  rate  proportional 
to  the  magnitude  of  both  of  these  quantities.  This 
heat  is  added  to  the  electron  species  directly  and 
then  passed  on  to  the  rest  of  the  gas  via 
collisional  processes.  The  reason  is  that  electrons 
carry  most  of  the  conduction  current  due  to  their 
low  mass.  In  a  steady  state  situation,  a 
compressible  gas  undergoing  Joule  heating  has 
an  elevated  electron  temperature  given  by  the 
following  relation  (Kerrebrock  [4]): 

!lZlL  =  JL(^KyM2o)1r1 

Tg  3  S' 

where  Te  and  Tg  are  the  electron  and  gas 
temperatures  respectively.  The  quantity  6’ 
denotes  a  mean  fractional  energy  loss,  and  is 
about  44  at  6000  K  in  air.  M  is  the  flow  Mach 
number,  cot  denotes  the  hall  parameter.  K  = 
E/(uB)  is  the  MHD  interaction  parameter,  and  y 
is  the  ratio  of  specific  heats  of  the  working 
medium. 

3.  Electron  beams:  Electron  beams  produce 
directed  ionization  in  gases.  They  provide 
convenient  means  to  generate  locally  high 
electrical  conductivity  at  relatively  low  flow 
temperatures.  The  electrical  conductivity 
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generated  by  an  electron  beam  is  usually 
obtained  from  the  following  expression: 


'  e1  ' 

i  (j>pre>' ] 

mnkc 

\  W,kjr  l  «  J 

where  e  is  the  electron  charge,  m  is  the  electron 
mass,  n  is  the  neutral  species  number  density,  kc 
is  an  electron  scattering  rate  constant,  kdr  (=  1.5e- 
7  cmA3/s)  is  a  dissociative  recombination  rate 
constant,  w*  is  the  ionization  constant  ( =  34  eV), 
jb  (about  0.1  A/cmA2)  is  the  current  density  in 
.  the  electron  beam,  Eb  is  the  energy  of  the 
electrons,  p  is  the  flow  density  and  Y  is  the 
electron  stopping  power.  The  energy  of  the  e- 
beams  can  be  varied  in  the  range  of  10  to  200 
keV,  to  provide  high  enough  electrical 
conductivity  at  moderate  flow  static  pressures. 
The  efficiency  of  electron  beams  diminishes  with 
increasing  pressure  due  to  the  inherent  nature  of 
the  electrons  to  recombine  with  greater  ease  in 
such  situations. 

4.  Seeding  with  easily  ionizable  alkali  salt: 
Traditionally,  MHD  power  generation  concepts 
have  involved  seeding  combustion  effluents  with 
some  easily  ionizable  salt,  such  as  Potassium 
Carbonate,  Cesium,  Sodium-Potassium 
solutions,  and  so  forth.  There  are  many  practical 
issues  involved  in  the  seeding  process,  involving 
multi-phase  flow,  deposition  of  seed  particles, 
conglomeration  of  seed  particles  near  cold  walls, 
and  the  resulting  electrode  shorting  and  so  forth. 
Nevertheless,  well-seeded  flows  exhibit  high 
electrical  conductivity  at  a  wide  range  of 
pressures  and  temperatures.  The  seed  material 
ionizes  essentially  in  an  equilibrium  manner.  The 
electrons  thus  produced  can  still  be  thermally  in 
nonequilibrium  with  the  rest  of  the  flow. 
Predicted  values  of  electrical  conductivity  in  the 
region  of  100  -  200  mho/m  can  be  derived  at 
temperatures  of  3000  K  and  flow  pressures 
between  1  and  10  atm. 

Vibrational-electronic  nonequilibrium  in  MHD 

Numerical  illustrations  of  supersonic  flow  in 
thermochemical  nonequilibrium  are  presented  here 
with  and  without  MHD  effects.  The  thermal  model 
used  here  is  the  Park  two-temperature  model  (Ref. 
[11]).  In  this  model,  the  vibrational-electron- 
electronic  energies  are  described  by  one  temperature, 
referred  to  as  Tv.  The  translational -rotational 
temperatures  are  assigned  the  symbol  T.  The  formal 
justification  for  this  model  may  be  found  in  Ref.[l  1]. 
This  model  is  useful  in  MHD  studies  of  air,  since  it 
provides  a  simple  framework  in  which  various  modes 


of  thermal  nonequilibrium  can  be  considered  in  a 
complex  mixture  of  gases.  The  strong  coupling 
between  the  electron  and  vibrational  modes  is  valid 
in  many  types  of  MHD  flows  of  importance.  It  is, 
however,  invalid  when  detailed  studies  of  the  anode 
sheath  is  undertaken. 

The  first  case  studied,  is  that  of  a  supersonic  flow  of 
air  (3000  m/s  at  T  =  3,000  K,  p  =  1  atm.).  The 
geometry  used  is  a  ramp  of  slope  lA  over  a  length  of 
0.8  m.  An  11  species  model  of  air,  including  N,  O, 
NO,  N2,  02,  N+,  0+,  NO+,  N2+,  02+  and  e-  is  used 
to  model  the  dissociation  phenomena.  Specific  heats 
are  deduced  from  curve  fits.  Vibrational  relaxation  is 
assumed  to  obey  the  Landau-Teller  theory.  Cases  in 
which  the  inflow  temperatures  T  and  Tv  are  both 
3,000  K  are  compared  with  the  situations  when  Tv  is 
forced  to  be  much  higher  at  the  inflow.  This  mimics  a 
situation  in  which  energy  is  deposited  in  the 
electronic  mode  by  means  of  some  suitable 
mechanism  at  the  inflow  of  a  supersonic  inlet.  Fig. 
[2]  shows  a  plot  of  the  thermal  relaxation  process 
across  a  leading  shock  wave.  It  is  seen  that  rapid 
energy  transfer  takes  place,  bringing  the  two 
temperatures  together  in  relatively  short  distance.  At 
higher  levels  of  Tv,  the  relaxation  process  is  faster 
due  to  the  increased  frequency  of  collisions.  The 
post-shock  increase  in  the  vibrational  energy  lags  the 
translational  energy.  This  is  due  to  the  fact  that  the 
vibrational  mode  experiences  the  shock  discontinuity 
only  through  relaxation  with  the  translational  mode. 

Fig.  [3]  shows  an  interesting  aspect  of  these  flows.  In 
the  nonequilibrium  situation,  the  electron 
concentration  has  been  artificially  set  to  3  orders  of 
magnitude  higher  than  the  equilibrium  value  in  this 
case.  Clearly,  the  electron  species  does  not  respond 
very  much  to  the  presence  of  a  shock  wave.  This  is 
evident  from  the  relatively  small  variation  of  the 
electron  density  across  the  shock  wave,  while  the 
Nitric  Oxide  and  other  species  exhibit  much  higher 
relative  variations.  This  indicates  that  the  electron 
species  responds  rather  sluggishly  to  shock  and 
contact  discontinuities.  In  flows  involving  MHD 
based  inlet  design,  this  manifests  as  a  loss 
mechanism,  since  the  electron  temperature  will  rise 
slowly  across  a  shock  wave,  thereby  making  the 
conduction  mechanisms  also  slower. 

The  second  case  studied  is  that  of  a  channel  with 
electric  and  magnetic  fields  applied  to  high 
temperature  supersonic  flow  of  air.  Air  flows  into  a 
constant  width  channel  at  3000  m/s,  6000  K,  1  atm., 
and  an  electric  field  of  strength  400  V/inch  is  applied 
along  the  y-axis  using  infinitely  segmented 
electrodes.  A  Hall  field  of  strength  about  1200  V/m  is 
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generated  in  the  x-direction  as  a  consequence. 
Constant  magnetic  field  of  strength  2.7  T  is  applied 
in  the  z-direction.  Fig.  [4]  shows  the  variation  of  the 
temperatures  T  and  Tv  with  and  without  the  MHD 
fields.  First,  the  lower  two  curves  show  the 
equilibrium  behavior  of  T  and  Tv  in  the  absence  of 
MHD.  There  is  a  slight  dip  in  their  value  at  the  inlet 
due  to  chemical  relaxation  of  gas  composition  to  the 
correct  equilibrium  value.  When  the  MHD  is  turned 
on,  both  temperatures  rise  almost  linearly  across  the 
length  of  the  channel.  They  also  maintain  a  uniform 
distance  between  them,  which  can  be  obtained  from 
the  closed  form  equation  almost  exactly.  This  serves 
as  a  validation  of  the  CFD  simulations.  Fig.  [4] 
shows  the  electron  and  NO+  mole  fractions.  It  is  seen 
that  with  or  without  MHD,  there  is  not  a  significant 
chemical  nonequilibrium  generated  at  this  pressure 
level.  The  electron  species  retains  charge  neutrality  in 
the  gas. 

Numerical  formulations  and  perspectives 

Maxwell  -Navier-Stokes  coupled  equation  set 
The  conservation  equations  governing  the  present 
applications  may  be  written  in  the  following  form: 

— + V  •  (Fx>  Fy ,  Fz)  +V  *  (Hx,  Hy,  Hz)  =  S 

dl 

in  which  the  fluxes  F  are  inviscid  or  convective, 
while  H  are  viscous  or  dissipative.  The  source  term  S 
comprises  of  the  production  of,  and  interaction 
between  the  various  flow  quantities.  This  equation 
may  be  expressed  in  the  integral  conservation  form  in 
a  control  volume  Q  bounded  by  dQ ,  as  follows: 

d  0+  (jfo  ds  +  H„  ds]=  js  dQ. 
cn  n 


The  fluxes  Fn  and  Hn  now  corresponding  to  normal 
face  fluxes.  The  various  components  in  this  equation 
are  described  below.  An  index  s  is  used  for  the 
various  gaseous  species  that  are  considered.  The  gas 
model  is  thought  to  comprise  of  the  components:  N, 
O,  K,  N2,  02,  NO,  A r,  0\  NO\  N2\  02\  tC,  e  .  In 
this  work  we  wish  to  consider  a  two  temperature 
model  of  the  flow.  The  reason  is  that  energy  addition 
by  MHD  at  high  temperatures  can  bring  about  a 
departure  from  equilibrium  between  the  various 
energy  modes.  In  cases  where  the  gas  is  weakly 
ionized  and  is  in  equilibrium,  small  perturbations 
from  the  equilibrium  state  may  be  modeled  based 
upon  the  Saha  equation  for  electron  number  density. 
In  the  following  equations,  we  will  use  the  two 
temperature  model  of  Park[l  1].  There  are  two  energy 


equations,  one  for  the  total  energy  conservation  and 
the  other  describing  the  vibrational-electronic  energy. 
The  following  matrices  describe  the  conservation 
laws: 
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The  transfer  of  energy  across  the  various  energy 
levels  is  contained  in  the  nonequilibrium  energy 
production  rate&>v.  This  term  contains  the  energy 
transfer  to  the  electronic  and  vibrational  modes  via 
collisions,  electron  pressure  gradient,  Landau-Teller 
type  of  relaxation  of  these  modes  with  the 
translational  and  rotational  modes,  and  release  in 
vibrational  and  electronic  energy  caused  by 
ionization  or  the  change  in  the  number  densities  of 
polyatomic  molecules.  The  expressions  used  to 
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model  these  features  may  be  found  in  various 
references  (e.g.:  Munipalli  [9,10]).  Radiative  energy 
addition  by  means  of  microwave  and  other  beams  can 
possibly  be  used  to  cause  an  increase  in  the  electron 
energy  and  ionization.  This  is  modeled  by  the  terms 
Qrad.  Energy  is  distributed  among  various  modes 
through  transport  processes.  Several  additional 
relations  are  required  to  complete  the  description  of 
nonequilibrium  MHD  gas  flow.  The  calculation  of 
physical  properties  is  described  in  the  next  sub¬ 
section.  Fluxes  F„  and  Hn  are  written  at  the  cell  face 
in  each  control  volume.  The  normal  components  of 
the  velocity  and  the  Magnetic  field  at  the  cell  face 


U„=V -n,Bn=Bn 

The  shear  stress  vector  rMis  related  to  the  flow 
velocity  gradients  and  the  coefficients  of  viscosity. 
Using  Stokes’  assumption  and  the  tensor  notation,  we 
may  write  the  shear  stress  as  a  flux  vector  of  the 
following  form: 


3  u  3c, 


T*=W 


The  two  temperature  model  used  above  may  be 
formally  justified  in  several  ways.  It  is  known  that 
the  electro-magnetic  fields  add  energy  directly  to  the 
“electron  pool”  of  the  flow.  Therefore,  it  is  essential 
to  model  the  electron  energy  as  a  distinct  thermal 
scale,  since  it  transmits  energy  to  other  energy  scales 
at  a  finite  rate.  The  other  modes  of  energy  storage, 
such  as  Vibrational,  Rotational  and  Translational 
energies,  may  be  resolved  using  the  arguments 
presented  by  Park  [11].  Equilibriation  between 
rotational  and  translational  modes  are  most  rapid  at 
all  temperature  ranges.  It  is  also  true  that  these  modes 


Temperature  (K) 


are  fully  excited  at  moderate  temperatures.  The  bulk 
of  the  vibrational  energy  is  stored  in  the  diatomic 
nitrogen,  since  it  is  the  species  with  the  largest  mass 
fraction  over  an  extended  temperature  range. 
Equilibriation  between  the  vibrational  energy  of 
Nitrogen  gas  and  the  electron  translational  energy  has 
been  studied  by  Lee [5].  In  the  temperature  range  of 
5,000  to  10,000  K  these  two  energy  levels  have  rather 
moderate  relaxation  times  and  in  other  words, 
equilibriate  most  rapidly.  This  in  essence  justifies  the 
use  of  a  two  temperature  model  as  described, 
comprising  of  T  and  Tv,  as  an  alternative  to  a  more 
elaborate  multi-temperature  model.  This  reduces 
stiffness,  computational  time,  required  modeling  data 
and  still  retains  a  fair  level  of  accuracy  in  the 
physical  modeling.  The  exact  trade-offs  in  this 
assumption  will  need  to  be  examined  after  a  realistic 
flow  configuration  is  chosen. 

Electric  Fields: 

In  the  above  formulation,  Hall  effect  has  been 
neglected.  The  wall  conductivity  distribution  could 
have  significant  influence  on  the  flow  of  Hall 
currents.  This  effect  will  be  reintroduced  in  these 
governing  equations  if  the  need  is  felt  at  a  later  stage. 
For  the  present,  by  rewriting  the  Ohm’s  law  relation: 
j  -  a  (  E  +  V  x  B  ),  and  using  Maxwell’s  equation: 
Vxfl  -  jjj\  we  obtain  the  following  expression  for  the 
electric  field: 


This  relation  will  also  be  used  in  applying  boundary 
conditions  on  the  electric  field  of  a  conducting  wall. 
In  the  case  of  a  perfectly  conducting  wall,  j  only  has 
a  normal  component.  In  an  insulator,  j  has  no  normal 
component.  If  the  magnetic  field  remains  steady,  the 
electric  field  may  be  expressed  in  terms  of  a 
potential.  This  would  cause  large  simplifications  in 
the  problem,  and  essentially  decouple  flow  equations 
from  Maxwell  equations. 

Transport  Properties 

In  the  study  of  high  temperature  gases,  the  transport 
properties  such  as  viscosity,  diffusion  coefficient, 
thermal  and  electrical  conductivities  are  evaluated 
from  the  momentum  transfer  collision  cross  sections. 
These  cross  sections  have  been  tabulated  for  the 
collisions  between  various  species  in  various  texts 
(see  for  example,  [8]).  The  expressions  used  to  derive 
some  of  these  quantities  are  written  out  below.  A 
sample  chart  showing  the  collision  cross  sections  of 
some  important  gases  with  the  electron  species  is 
shown  in  Fig.[l]. 


Figure  1:  e‘-collision  cross  sections 
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Using  w4,to  denote  the  particle  mass  in  kg  of  a 
species  s,  and  A(4J?and  A^  to  denote  collision 
integrals  between  the  species  s  and  r  and  the 
following  definition  for  the  molar  concentration: 


Ys  = 


Ps 

pMs 


Coefficient  of  viscosity  is  given  by: 


*=Z 
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Translational  thermal  conductivity  is  given  by  (using 
kio  denote  Boltzmann’s  constant): 


Boundary  Conditions: 

Boundary  and  initial  conditions  for  this  equation  set 
must  be  specified  accurately  in  order  to  prevent 
numerical  instabilities  and  non-physical  solution 
behavior.  The  flow  variables  such  as  velocity,  density 
and  pressure  can  be  specified  from  an  equilibrium  or 
perfect  gas  estimate  at  the  prescribed  flight  condition. 
Boundary  layer  values  of  species  densities, 
temperature  and  turbulence  energy  and  intensity  must 
be  given  as  initial  conditions.  For  supersonic  and 
super-alfVenic  flows,  information  propagates  only  in 
the  downstream  direction.  The  inflow  variables  are 
then  fixed  and  outflow  variables  are  copied  from 
within  the  domain.  For  subsonic  or  sub-alfVenic 
flows,  characteristic  relations  must  be  used  in  order 
to  propagate  information  either  into  or  out  of  the 
domain. 


n'  4  *5-  2^£>(r)+3.54n  A^>(re) 

r*e~ 

Rotational  and  electronic  conductivities  are  defined 

in  a  similar  manner,  using  A(1)  in  the  place  of  A(2) . 
The  vibrational  thermal  conductivity  is  assumed  to  be 
equal  to  the  rotational  conductivity.  The  species 
diffusion  coefficients  are  given  by: 

(Z  r$Ms{\-Msys) 

V  YrP^viT) 

h  W 

The  electrical  conductivity  of  a  partially  ionized  gas 
is  given  by  the  following  expression,  using  ne,  nt  and 
ns  to  denote  the  number  densities  of  the  electrons, 
ions  and  unionized  species  respectively: 


The  wall  (airframe  surface)  may  be  set  to  be  a  no-slip 
non-catalytic  isothermal  boundary.  This  prescribes  a 
constant  wall  temperature  distribution,  zero  normal 
derivatives  for  species  densities  and  zero  velocity  at 
the  wall.  Electromagnetic  boundary  conditions  may 
be  more  involved,  depending  on  the  choice  of  wall 
material.  For  a  perfect  conductor,  the  normal 
component  of  the  magnetic  field  is  zero  and  the  jump 
in  the  electric  field  has  no  tangential  component.  For 
a  perfect  insulator,  the  gradient  of  the  magnetic  field 
is  zero  at  the  surface. 

Space/time  integration:  Numerical  Scheme 
A  first  order  accurate  explicit  finite  volume  scheme 
to  solve  the  governing  equations  shown  earlier  may 
be  written  for  a  cell  bearing  the  index  i,  volume  O,, 
bounded  by  with  edge  lengths  As  as: 

ur=u; 

an, 


Z"sA?M>  +  3-9"; 
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log 
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In  the  above,  e  is  the  electron  charge  and  ce  is  the 
mean  thermal  velocity  of  the  electrons.  The  Hall 
parameter  is  now  be  deduced  from: 


The  fluxes  F  and  H  above  must  be  defined  at  the  cell 
faces  in  an  appropriate  manner  in  order  to  be  able  to 
capture  shock  waves  accurately.  Though  the  present 
set  of  equations  has  not  been  solved  using  such  a 
scheme  in  any  of  these  references,  each  reference 
validates  a  portion  of  the  complete  set  of  equations 
presented  here.  It  is  anticipated  that  strong  shock 
waves  near  blunt  body  regions,  as  well  as 
predominantly  viscous  regions  of  a  flow  can  cause 
numerical  difficulties  with  this  implementation  of 
Roe’s  scheme.  In  order  to  overcome  these 
difficulties,  other  flux  functions  will  be  built  into  this 
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solver.  Kinetic  flux  vector  splitting  schemes  have 
been  successfully  tested  on  ideal  fluid  MHD 
problems  by  Xu[14],  Linde  [6].  The  inclusion  of 
kinetic  schemes  is  computationally  simpler.  A  full 
description  may  be  found  in  [6].  However,  these 
schemes  are  very  diffusive  and  less  accurate  in 
comparison  with  flux  difference  split  strategies, 
particularly  for  shear  and  contact  waves.  A 
combination  of  both,  as  described  in  [6]  could  also 
yield  a  viable  solution  procedure. 

Due  to  inherent  limitations,  subsonic  flow  features 
are  heavily  damped  by  upwind  schemes.  This  causes 
degradation  in  viscous  flow  solutions.  Wall  shear, 
heat  transfer  rates  and  turbulent  flow  features  are 
affected  by  this  dissipation.  Higher  order  corrections 
must  therefore  be  made  in  order  to  resolve  this  issue. 
An  alternate  numerical  procedure  to  the  above,  is 
described  by  Agarwal  and  Augustinus  [1].  Here,  a 
Runge-Kutta  time  stepping  finite  difference  scheme 
with  TVD  corrections  has  been  used.  A  second  order 
symmetric  TVD  model  devised  by  Davis  and  Yee  has 
been  proposed.  The  advantages  of  this  approach  over 
Roe-type  solvers  for  stiff  equations  with  coupled 
flow  physics  is  yet  unknown.  During  the  course  of 
this  study,  by  controlling  the  effects  of  the  variables 
such  as  Magnetic  field  intensity,  electrical 
conductivity  and  so  forth,  it  is  possible  to  compare 
and  converge  upon  a  numerical  scheme  that  will 
address  all  of  the  concerns  described  above. 


The  KFVS  scheme  has  been  applied  to  the  1-D  shock 
tube  with  finite,  zero  and  infinite  electrical 
conductivity.  For  the  results  shown  here,  a  modified 
Brio  and  Wu  [2]  problem  has  been  studied,  where  the 
driver  and  driven  conditions  are: 


Pressure 

l.e6  l.e5 

(Pa) 

Density 

10.  1.25 

(kg/mA3) 

Bx 

0.75  0.75 

(T) 

By 

1.  -1. 

(T) 

T 

2. 

cr 

0.,  1000.,  oo 

(fim)-1 

At 

4.e-5 

(s) 

At  /  Ax 

4.e-4 

(s/m) 

400  mesh  points  were  used  for  the  shock  tube  and  the 
solution  was  plotted  after  200  time  steps.  The 
parameter  rj  from  Ref.  [14]  was  fixed  at  0.5  for  all 
the  cases.  A  CFL  number  of  about  0.8  results  from 
the  above  selection.  Fig. [6]  shows  a  plot  of  pressure, 
velocity  and  Magnetic  field  for  the  various  situations 
of  interest. 

It  has  been  observed  that  a  close  correlation  exists 
between  the  magnitude  of  the  electrical  conductivity 
and  the  permissible  value  of  Ax.  This  is  due  to  the 
presence  of  a  second  derivative  with  respect  to  x  on 
the  right  hand  side  of  the  governing  equations  when 
or  is  nonzero.  The  extra  Ax  term  violates  scaling  with 
respect  to  the  conventional  CFL  number  and  presents 
itself  as  a  singular  perturbation  parameter. 


Results 

Work  performed  on  decoupled  Navier-Stokes  - 
Maxwell  equations  for  high  temperature  seeded  and 
unseeded  air  has  been  reported  earlier  (Munipalli  et 
al.  [9,10]).  In  this  paper,  we  present  simple 
applications  of  the  coupled  equation  set,  the  nature  of 
this  coupling  and  some  critical  issues  related  to 
boundary  conditions  and  numerical  stability.  This 
paper  represents  an  initial  step  towards  the  attainment 
of  a  large  scale,  unstructured  grid  parallel  code 
capability  to  simulate  MHD  problems. 

Three  schemes  are  examined  here.  The  first  is  a 
Kinetic  Flux  Vector  split  scheme  as  presented  by  Xu 
[14]  for  one-dimensional  unsteady  cases.  The  second 
is  the  Jameson-Schmidt-Turkel  [3]  central  difference 
based  Runge-Kutta  time-stepping  scheme,  used  in 
two-dimensional  problems.  Roe’s  scheme,  as 
presented  by  Linde  [6],  is  the  third  scheme  used. 
Finite  volume  discretizations  and  dimensional  flow 
quantities  are  used  thoughout. 

One-Dimensional  Flows 


The  Neumann  stability  of  this  system,  therefore 
follows  that  of  the  viscous  Burgers’  equation,  and 
two  parameters  must  be  tracked  independently: 


Where  u+cf  is  the  fastest  wave  speed  at  a  cell 
interface.  The  stability  criterion  is: 


v2  <2 r 


Results  presented  in  Fig.  [6]  show  a  conductivity  of 
1000  (fim)-1,  which  was  arbitrarily  chosen.  The 


scheme  however,  functions  well  for  this  case. 


Nature  of  the  NS -Maxwell  coupling 
Lorentz  force  and  Joule  heating  terms  may  be  added 
trivially  to  most  existing  Euler  or  Navier  stokes  codes 
in  situations  where  these  terms  are  small  and  the 


induced  magnetic  field  is  negligible.  While  this 
assumption  is  valid  in  several  cases  of  engineering 
importance,  a  truly  coupled  solver  is  important  when 
electrical  conductivity  increases  and  the  length  scales 
available  for  MHD  interaction  are  large.  In  order  to 
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assess  the  relative  importance  of  coupling  vs. 
decoupling  the  Navier-Stokes  and  Maxwell 
equations,  the  Brio-Wu  test  case  described  above  has 
been  simulated  with  various  models. 

Firstly,  an  Euler  code  based  on  Roe’s  scheme  (made 
second  order  accurate  using  MUSCL  type  of 
upwinding,)  was  modified  to  include  MHD  terms. 
The  simplest  manner  in  which  the  one  dimensionality 
of  the  problem  may  be  non-trivially  retained  is  to 
assume  that  the  velocity  is  purely  uni-directional  (x) 
and  the  magnetic  field  is  applied  in  a  perpendicular 
direction  (z).  In  this  case,  as  shown  in  Fig.[7],  the 
source  terms  dissipate  the  velocity  profiles  and  damp 
out  the  discontinuities  in  the  pressure  profile.  When 
the  electrical  conductivity  is  finite,  this  is  indeed  the 
expected  effect. 

A  similar  experiment  was  conducted  with  the  KFVS 
based  code  described  earlier.  Magnetic  field,  now,  is 
oriented  in  the  x-y  plane,  and  is  kept  constant  through 
out  the  channel  as  Bx  =  0.75  T  and  By  =  1  T.  Cases 
were  run  with  and  without  enabling  the  update  of 
magnetic  field  components  in  the  code.  As  can  be 
seen  from  Fig.  [8],  significant  departures  (upto  100% 
and  more!)  may  be  observed  even  in  cases  with  a 
finite  conductivity.  The  case  where  B  was  held 
constant  is  seen  to  retain  a  wave  structure  similar  to 
perfect  gas,  while  having  a  distinctly  increased 
distance  between  the  leading  shock  wave  and  the 
contact  surface.  This  aspect  will  be  the  subject  of  a 
future  paper  and  was  discussed  briefly  in  [9], [10]. 

Two  Dimensional  Flows 

In  the  2-D  cases,  both  Roe’s  approximate  Riemann 
solver,  as  well  as  the  Jameson-Schmidt-Turkel  [3] 
strategy  have  been  adopted.  The  latter,  primarily  due 
to  its  simplicity  in  coding,  and  available  experience 
with  its  use  in  fluid  mechanics  [3],  electromagnetics 
(Shankar  et  al.  [11]),  and  recently,  in 
magnetohydrodynamics  (Agarwal  et  al  [1]).  Second 
order  damping  was  used  in  supersonic  flow 
problems,  and  the  MHD  pressure: 

*  B2 

p  -  p  +  — 

F  F  2  ju 

has  been  used  to  compute  the  damping  coefficient. 
Normal  gradients  of  this  pressure  has  been  set  to  zero 
at  the  wall,  and  the  magnetic  field  component 
tangential  to  the  wall  has  been  fixed  as  zero. 

Preliminary  results  for  the  case  of  a  supersonic  flow 
past  a  ramp  are  presented  here.  Thin  layer  Navier 
Stokes  -  Maxwell  equations  have  been  solved  for  air 
(y  =  1.4,  R  =  277  J/(mol.K)),  velocity  derivatives 
along  the  wall  are  assumed  to  not  contribute  to  the 


shear  stress  terms.  No  special  treatment  for  making 
the  magnetic  field  divergence  free  has  been  adopted. 
Electrical  conductivity  is  assumed  to  be  infinite  for 
the  present. 

In  the  Roe  solver,  however,  the  source  term  proposed 
by  Powell  [12]  has  been  used  to  overcome  the  crisis 
created  by  the  wave  structure.  An  eight  wave  system 
is  used  in  the  manner  presented  by  Linde[6].  Central 
differences  have  been  used  to  compute  Roe 
averaging.  (A  more  complex  density  based  averaging 
was  also  used,  but  the  difference  was  unnoticeable 
for  the  cases  studied.)  Firstly,  the  Brio-Wu  shock 
tube  problem  was  simulated  using  a  two-dimensional 
shock  tube  with  periodic  wall  conditions.  Results, 
shown  in  Fig.  [9]  match  closely  with  1-D 
simulations.  Magnetic  field  in  the  x-direction,  which 
is  fixed  at  0.75  in  the  1-D  7  wave  system,  is  seen  to 
vary  about  3%  depending  on  the  type  of  wall 
boundary  conditions  used.  This  fluctuation  is 
sensitive  to  the  Powell-type  source  terms  and  can 
change  significantly  if  they  are  absent. 

Flow  past  a  5  degree  ramp  at  Mach  6  with  applied 
magnetic  field  in  the  y-direction  was  studied.  The 
presence  of  fast  and  slow  magneto-acoustic  waves 
are  easily  seen  in  these  results.  Fig.  [10]  shows  the 
induced  field  vectors.  Near  the  body  surface,  the 
induced  field  acts  in  such  a  way  as  to  annul  the 
tangential  magnetic  field  component.  The  normal 
component  of  the  magnetic  induction  vector  H  is  held 
continuous  across  a  boundary.  Magnetic  field  B  is 
computed  from  the  relative  values  of  permeability  in 
the  two  media,  which  was  taken  as  14  in  the  present 
case.  The  formation  of  magnetosonic  shocks  for 
various  values  of  applied  field  is  shown  in  Fig.  [11]. 
The  slow  shock  increases  in  strength  as  B  increases 
and  is  capable  of  overtaking  the  fast  shock 
eventually.  A  comparison  with  the  Jameson-type 
solver,  presented  in  Fig  [12]  shows  minor  oscillations 
which  may  be  controlled  with  the  appropriate 
modification  of  the  dissipation  terms  (which  are 
currently  pressure  dependent.)  Perhaps  a  eigen  vector 
based  limiter  would  be  more  appropriate. 

Code  Development  Environment 
An  existing  code  environment  for  Computational 
Electromagnetics  is  being  modified  to  incorporate  the 
algorithms  proposed  here.  The  code  environment  is 
capable  of  running  in  a  parallel  mode  using  MPI  and 
is  being  used  extensively  on  PC  clusters  running 
LINUX. 

Conclusion 
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The  current  thinking  towards  the  development  of 
a  comprehensive  real  gas  MHD  solver  on 
unstructured  grids  with  thermochemical 
nonequilibrium  has  been  described  in  this  paper. 
A  two  temperature  nonequilibrium  decoupled 
MHD  solver  was  developed  in  a  prior  work, 
which  will  be  interwoven  into  a  Navier- Stokes  - 
Maxwell  coupled  solver.  Schemes  based  on 
approximate  Riemann  solvers,  as  well  as 
Jameson-type  schemes  have  been  seen  to  be 
effective  in  the  solution  to  the  coupled  equation 
set.  Some  interest  also  exists  in  applying  the 
more  recently  developed  Kinetic  Flux  Vector 
Splitting  scheme  to  the  coupled  equation  set. 
Applications  to  hypersonic  inlet  mass  capture  and 
shear  layer  stability  (Kelvin-Helmholtz  type)  are 
being  studied  currently  using  these  models. 
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Figure  2  :  Thermal  relaxation  process  at  T  =  3000 
K  and  Tv  =  6000  K 


Figure  3:  e-  and  NO  distributions  for  ramp 
flow  with  and  without  thermo-chemical 
equilibrium 
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Figure  4  :  Nonequilibrium  Temperature  variation 
along  an  unseeded  air  MHD  accelerator  with  and 
without  the  MHD  fields 
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Figure  5:  Variation  of  electron  and  NO+  mole 
fractions  along  a  nonequilibrium  unseeded  air 
MHD  accelerator 


Figure  6:  Flow  quantities  for  the  modified  Brio  and  Wu  problem 
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Figure  7:  MHD  terms  included  as  sources  in  an  Euler  equation  solver 
( Brio-Wu  shock  tube  problem  ) 


Figure  8:  The  effect  of  coupling  the  NS-Maxwell  equations  for  moderate  induced  fields 
( Modified  Brio-Wu  shock  tube  problem  ) 
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Figure  11  :  Formation  of 
fast  and  slow  magnetosonic 
shocks  in  2-D  ramp  flow 
problem 
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Figure  12:  Comparison  between  Roe’s 
scheme  and  Jameson-type  central 
difference  scheme  for  an  MHD  oblique 
shock  wave  problem 
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Abstract 

In  this  paper  we  describe  the  development  strategy  for  a 
comprehensive  unstructured  grid  based  simulation 
system  for  magnetohydrodynamics  (MHD)  flow  over 
large  complex  geometries.  Earlier  work  has  reported  the 
successful  use  of  modern  CFD  techniques  in  simulating 
MHD  based  devices  ([14], [15]).  Here,  we  focus  on 
hypersonic  inlet  design  augmented  with  the  use  of 
MHD.  Boundary  layer  effects,  inlet  mass  capture,  and 
other  flow  nonuniformities  are  discussed.  A  hierarchical 
treatment  of  physics  in  MHD  solvers  described  here, 
provides  means  to  both  validate  numerics  and  physics  in 
a  step-by-step  manner,  as  well  as  develop  confidence 
estimates  for  simpler  models.  Real  gas  effects  include 
finite  conductivity,  thermal  nonequilibrium,  charge 
separation,  chemical  reactions,  and  other  diffusive  and 
dissipative  effects. 

Hierarchical  MHD  flow  physics 

Limitations  of  Source  term  approach:  The  most 
straightforward  method  in  which  MHD  effects  may  be 
introduced  into  an  existing  CFD  code  is  by  the  addition 
of  source  terms  corresponding  to  the  Lorentz  force  in 
the  momentum  equation,  and  Joule  heating  in  the 
energy  equation  with  an  initially  prescribed  distribution 
of  electric  and  magnetic  field.  Such  an  approach  is 
simple  to  implement  in  an  existing  fluid  flow  code.  The 
applicability  of  this  method  is  limited  by  the  size  of  the 
magnetic  field  that  is  induced  in  the  flow  by  virtue  of 
the  currents  generated.  For  large  values  of  induced 
fields,  the  magnetic  field  must  be  computed 
simultaneously  at  each  time  step.  Pressure  tends  to  be 
underpredicted  by  such  a  source  term  based 
implementation. 

Constant  B.  no  electric  field  present:  For  a  three 
dimensional  velocity  field  in  which  the  magnetic  field 
has  all  the  three  components,  Ohm’s  law  neglecting  the 
Hall  effect  and  electric  field,  may  be  written  as: 


J  =  a(fx  B)-a{{vBz  -wBy) i  + 

(i wBx-uBz  ) j  +  (uBy  -  vBx  )k  }  [  1  ] 

where  u  ,v  and  w  denote  velocity  components.  This 
relation  may  be  written  for  simplified  2-D  situations 
where  the  magnetic  field  may  have  components  in  some 
preferential  directions  alone.  This  is  perhaps  the 
simplest  manner  in  which  to  include  the  effects  of  the 
magnetic  field  and  also  the  least  accurate.  For  large 
magnetic  Reynolds  numbers  as  well  as  for  cases  where 
electric  current  and  Hall  effect  are  important,  this 
relation  cannot  be  used. 


Constant  B.  with  spatially  varying  electric  field:  The 
general  Ohm’s  law  including  an  electric  field,  is  given 

by  J  =  cr^  +  V  xB)  which,  for  situations  in  which  the 

electric  field  may  be  deduced  from  a  potential  (which 
always  exists  when  the  magnetic  field  is  constant  in 
time,)  mav  be  rewritten  as: 

J  =  ar(-V<p  +  VxB)  [2] 

In  order  to  determine  this  potential,  the  divergence  of 
the  above  equation  which  is  zero,  by  the  law  of 
conservation  of  charge,  is  taken,  yielding  the  following 
equation: 

V20>  =  V  •  (rx  fi)  [3] 


Various  simplifications  may  now  be  applied  to  this 
equation.  If  a  space  varying  magnetic  field  acting  in  the 
z-direction  is  applied  to  a  two-dimensional  velocity 
field, 


vV  = 


8  (vB. ) 
8x 


d(uB.) 

dy 


[4] 


When  the  electrical  potential  is  computed,  it  is  readily 
substituted  in  eq.  [2]  to  obtain  the  value  of  J  and  thus, 
the  Lorentz  force  JxB  may  be  deduced.  This  set  of 
equations  must  be  augmented  by  the  appropriate 
boundary  conditions  described  in  the  following  section. 
These  equations  may  be  written  out  trivially  for  a 
general  three  dimensional  situation  of  current  interest. 
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Induced  magnetic  field  in  one  direction:  The  situations 
when  the  applied  magnetic  field  is  oriented  along  a 
fixed  direction,  and  the  induced  field  is  dominant  in  one 
direction  have  been  considered  by  several  authors  (e.g., 
Morley  [12],  Shishko  [18]).  Consider  the  case  when  the 
induced  field  Bx  is  in  the  x-direction  and  the  applied 
field  (0,By,Bz)  is  oriented  in  the  y-z  plane.  The 
momentum  equation  may  be  rewritten  in  terms  of  the 
induced  field,  as  opposed  to  the  current  as  used  earlier. 
This  and  the  time  advancement  of  the  induced  field,  are 
governed  by: 


dpu  d(/w2  +  P-RbBllfim) 


8t 


Ry 


R, 


V 


dBx 

dy 


+  B, 


8x 

gg, 

z  8z 

J 

du 


\ 


y8y  1  8z 


8t 


V2w 

Re 

1 

Re*  % 


[5] 


-V2fl 


This  equation  without  the  time  derivatives  and  assuming 
a  fully  developed  channel  flow  in  the  x-direction  has 
been  applied  with  insightful  results  in  [12], [18].  Present 
work  provides  an  opportunity  to  study  time  varying 
induced  fields  for  general  cases.  The  formulation  for  a 
general  3-D  situation  where  all  the  components  of  the 
induced  field  exist,  (and  the  definition  of  the  non- 
dimensional  quantities  used  above),  has  been  presented 
in  Ref.  [14]. 


Integrated  NS-Maxwell  equations:  A  full  complement 
of  Navier-Stokes  and  Maxwell  equation  set  has  been 
solved  at  HyPerComp  for  problems  of  interest  in 
compressible  flow  (Ref.  [14]).  The  principal  difficulty 
in  solving  this  set  of  equations  is  to  maintain  a 
divergence  free  magnetic  field  (Brackbill  and  Barnes 
[1]).  Characteristics  based  techniques  exist  in 
compressible  flow  solvers,  designed  to  avoid  this 
problem.  In  incompressible  flows,  one  must  post¬ 
process  the  magnetic  field  to  remove  the  effects  of  non¬ 
zero  divergence.  This  typically  involves  the  solution  to 
a  Poisson-type  equation  and  may  be  computationally 
expensive.  However  the  equation  using  the  electric 
potential  also  has  this  difficulty.  The  Maxwell  equations 
are  written  in  terms  of  the  magnetic  field  components. 


High  temperature  effects:  Hypersonic  MHD  is  virtually 
characterized  by  an  abundance  of  thermochemical 
nonequilibrium  processes.  High  temperatures 
encountered  in  such  flows  cause  significant  variations  in 
the  gas  specific  heats,  chemical  composition,  viscosity, 
thermal  and  electrical  conductivity  and  the  functional 
dependence  of  the  speed  of  sound  on  flow  properties. 
Typical  equilibrium  and  nonequilibrium  models  have 
been  presented  by  Munipalli  [13],  Munipalli  et  al.  [15] 
for  MHD  accelerators.  Joule  heating  transfers  energy 


preferentially  to  the  free  electron  translational  modes. 
By  virtue  of  collisions,  this  energy  is  imparted,  in  order 
of  ease,  to  the  internal  electron  modes,  the  vibrational 
modes  (if  present),  rotational  modes  and  translational 
modes.  While  a  detailed  balancing  of  energy  is  needed 
in  high  temperature  mixture  of  gases  (such  as  air),  we 
have  in  the  past  experimented  with  the  Park  two- 
temperature  model  with  fair  comparisons  with  expected 
semi-analytical  results. 

Nonequilibrium  treatment  near  walls:  The  effects  of 
nonequilibrium  in  MHD  are  most  pronounced  near  the 
electrode  region.  In  the  core  of  the  flow,  vibrational  and 
electron  modes  are  tightly  coupled,  and  a  two 
temperature  model  proposed  by  Park  [16],  as  evidenced 
form  the  work  of  Lee  [8]  may  be  adequate.  .Near  an 
electrode  wall,  the  vibrationally  temperature1  quickly 
cools  to  the  boundary  value,  while  the  electron  mode 
loses  temperature  much  more  gradually,  well  into  the 
sheath  region,  in  the  vicinity  of  a  few  mean  free  paths 
from  the  wall.  This  brings  up  a  need  for  a  three 
temperature  model  in  order  to  resolve  the  thermal 
gradients  near  an  electrode.  Three  temperature  models 
such  as  those  used  by  Sleziona  et  al.  [19]  are  being 
currently  pursued. 

Numerical  Difficulties 

Problem  at  the  research  level  abound  in  both  physical  as 
well  as  numerical  understanding  of  MHD  processes. 
Some  of  these  are  described  here. 

Divergence  free  magnetic  fields:  Much  interest  in  recent 
times  has  led  to  the  evolution  of  schemes  which 
preserve  the  divergence  free  nature  of  the  magnetic 
field.  While  a  non-zero  divergence  of  magnetic  field  is  a 
serious  error,  it  is  unclear  from  our  current  experience, 
if  this  problem  has  as  much  a  bearing  in  the  simulation 
of  hypersonic  flows,  where  the  electrical  conductivity  is 
rather  low  as  compared  to  plasmas  encountered  in 
nuclear  physics  or  astrophysics. 

Accurate  modeling  of  electron  temperature:  While 
numerical  stiffness  dominates  the  simulation  of  finite 
rate  chemistry  processes  and  thermal  relaxation,  the 
near-electrode  region  presents  a  numerical  challenge  of 
a  different  nature.  The  length  scales  in  which  important 
transport  processes  occur  near  an  electrode  are  several 
orders  of  magnitude  smaller  than  all  other  scales  present 
in  the  flow.  In  ref.  [3]  a  solution  matching  procedure  has 
been  developed,  by  which  the  wall  processes  may  be 
analytically  reconstructed.  This  has  been  successfully 
applied  to  the  case  of  an  arc-heater,  and  it  is  believed 
that  a  similar  approach  may  be  brought  to  bear  in 
hypersonic  MHD. 
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Hartmann  number  limitations:  Large  Hartmann 

numbers  are  encountered  routinely  in  the  flows  of  liquid 
metals  in  nuclear  fusion  and  many  modem  metal  casting 
processes.  In  such  situations,  flow  velocity  gradients 
near  the  wall  become  several  times  enhanced,  and 
frequently  lead  to  numerical  instability.  This  problem 
remains  intractable  to  a  large  extent  to  date.  Leboucher 
[7]  presents  an  approach  by  which  wall  functions  and  a 
staggered  mesh  treatment  may  alleviate  this  concern  in 
incompressible  flows.  Fortunately,  this  does  not  seem  to 
be  a  major  issue  in  aerodynamics  due  to  the  low  values 
of  electrical  conductivity. 

Treatment  of  wall  boundary  conditions 


Fluid  Quantities:  The  wall  (airframe  surface)  may  be 
set  to  be  a  no-slip  non-catalytic  isothermal  boundary. 
This  prescribes  a  constant  wall  temperature  distribution, 
zero  normal  derivatives  for  species  densities  and  zero 
velocity  at  the  wall. 


Electrode  surface:  At  high  voltages  in  weakly  ionized 
gases,  a  sheath  layer  exists  close  to  electrode  surfaces. 
This  is  a  free-molecular  region  in  which  charge 
separation  can  exist  and  accounts  for  a  potential  drop,  as 
well  as  two-dimensional  nonuniformities  in  the  current 
distribution.  The  length  scale  of  the  sheath  layer  is  often 
several  orders  of  magnitude  smaller  than  the  boundary 
layer  thickness,  and  therefore,  one  must  take  recourse  to 
analytical  corrections  for  the  electromagnetic  effects 
occurring  in  the  sheath.  The  sheath  comprises  of 
electron-ion  recombination  processes  which  cease  when 
no  further  net  mass  influx  of  electrons  is  possible  due  to 
the  effective  potential  drop.  This  results  in 
(recombination)  energy  release,  which  is  applied  to  the 
energy  equation: 


q  Is  4 


8  kT 


nm. 


1/2 
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where  the  index  s  refers  to  ionized  particles,  and  1  refers 
to  the  ionization  energy.  Electron,  Ion  and  wall 
temperatures  are  identical  in  the  sheath.  The  difference 
in  the  mass  of  electrons  and  ions  entering  the  sheath 
layer  causes  a  potential  drop  given  by: 

kT . 


=  —  tog 
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These  processes  take  place  in  a  layer  whose  thickness  is 
of  the  order  of  the  Debye  length,  which  is  somewhat 
smaller  than  the  mean  free  path  of  electrons  in  the 
situations  being  considered.  Since  this  distance  is  too 
small  to  be  resolved  in  a  computational  mesh,  this 
relation  is  applied  as  a  boundary  condition. 


Effect  of  segmentation  of  electrodes:  In  applying 
current  boundary  conditions  on  a  wall  comprising  of 


segmented  electrodes,  certain  simplifications  must  be 
made.  If  the  computational  mesh  can  resolve  each 
electrode  sufficiently  (which  is  extremely  difficult  to 
accomplish  in  real  problems,)  one  would  set  the  normal 
component  of  current  to  zero  in  the  insulating  region 
and  the  tangential  gradient  of  the  electric  potential  to 
zero  on  a  conducting  surface.  However,  if  an  infinitely 
segmented  electrode  arrangement  is  assumed,  it  is 
possible  to  deduce  the  following  condition  on  the  wall 
(Munipalli  [29]): 

=  [8] 

dy  dx 

where  y  denotes  the  direction  normal  to  the  wall  and  x 
along  the  wall,  <p  is  the  electric  potential  and  (3  is  the 
Hall  parameter.  This  expression  seems  adequate  for 
practical  applications.  As  an  example  the  potential 
distribution  in  a  slightly  expanding  channel  from  AEDC 
experiments  in  the  1960s  has  been  computed  and  seen 
to  match  almost  exactly  with  test  data.  Results  are 
shown  in  Fig.[l].  117  electrode  pairs  are  installed  on 
parallel  walls  and  a  magnetic  field  is  applied  in  a 
perpendicular  direction.  While  the  potential  difference 
in  each  electrode  pair  is  fixed  by  connecting  a  battery 
across  it,  the  potential  difference  along  the  channel  is 
derived  by  using  this  condition  in  conjunction  with  a 
flow  solver. 

Magnetic  Field:  At  a  perfectly  conducting  wall,  normal 
component  of  the  magnetic  induction  H  is  continuous 
across  the  wall.  Frequently  in  computations,  the 
tangential  component  of  this  field  is  set  to  zero,  which  is 
an  approximation  for  slender  bodies.  The  magnetic  field 
strength  inside  the  solid  body  would  be  obtained  by 
scaling  the  gas  magnetic  field  strength  by  the  ratio  of 
the  magnetic  permeability  in  the  two  media.  The 
induced  magnetic  field  vectors  in  a  supersonic 
conducting  flow  past  a  ramp  is  shown  in  Fig.  [3].  Here, 
a  vertical  magnetic  field  is  applied,  which  is  nearly 
cancelled  tangential  to  the  wall  by  a  tangential  induced 
field  in  the  opposite  direction. 

Four  principal  methods  in  which  MHD  boundary 
conditions  are  being  implemented  in  our  work,  are 
explained  below. 

(a)  No  induced  field \  insulating  walls:  This  is  the 
simplest  situation  among  the  cases  of  interest.  No 
BC  is  needed  on  the  magnetic  field,  while  a 
Neumann  condition  is  applied  to  the  electric 
potential. 

(b)  Induced  field  with  insulating  walls:  Here,  the 
normal  component  of  the  induced  magnetic  field 
strength  vector  H  is  often  assumed  to  be  continuous 
across  a  solid  wall.  The  boundary  condition  on  the 
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normal  component  of  the  magnetic  field  is  written 
as  (where  p  denotes  the  magnetic  permeability): 

^ri, wall  _  Mm, wall 

Bn,  fluid  Mm,  fluid 

(c)  Thin  conducting  walls :  This  boundary  condition 
may  be  formulated  with  respect  to  either  an  electric 
potential,  if  there  is  no  induced  field,  or  an  induced 
magnetic  field,  if  present.  In  the  first  case,  when  the 
wall  is  conducting,  if  the  magnetic  field  is  assumed 
to  be  oriented  along  the  z-direction,  and  the  wall  in 
the  x-y  plane,  if  the  number  c  denotes  the  ratio  of 
the  electrical  conductivity  of  the  wall  to  that  of  the 
fluid  medium,  the  following  boundary  relation 
applies: 

j-k  =  -V<p-k  =  -Vxy-(ctVxy<p)  [10] 

where  t  denotes  the  non-dimensional  thickness  of 
the  wall.  The  gradient  operators  above  are 
evaluated  in  the  x-y  plane.  This  equation  provides 
for  the  continuity  for  the  flow  of  currents,  and 
provides  a  wall  potential  distribution  matching  the 
electric  field  in  the  flow.  A  similar  condition  may 
be  evolved  for  the  case  of  induced  fields,  and  for 
the  case  described  earlier  (induced  field  in  the  x- 
direction,  applied  field  in  the  y-z  plane,)  the 
following  relation  may  be  written  on  a  wall  placed 
in  the  x-y  plane: 

dBx  n  A 

Ct-TL+Bx=  0  [11] 

dz 

In  reality,  all  of  the  above  are  vector  relations,  and 
may  be  written  out  for  an  arbitrarily  oriented  mesh 
structure. 

The  formation  of  Hartmann  layers  perpendicular  to 
the  applied  magnetic  field  near  the  wall  regions 
distinguishes  MHD  flows  from  traditional  fluid 
boundary  layers.  Hartmann  layers  provide  a  path 
for  current  lines  from  the  core  flow  to  close.  The 
thickness  of  the  Hartmann  layer,  6,  is  set  equal  to 
the  reciprocal  of  the  Hartmann  number.  For  large 
Ha,  the  Hartmann  layers  are  so  thin  that  it  is 
difficult  to  resolve  them  numerically.  Leboucher 
[47]  reported  success  in  an  integral  treatment  of  the 
Hartmann  layers  and  obtained  stabilized  solutions 
for  Hartmann  numbers  up  to  1000.  This  involved 
adding  5  to  the  effective  thickness  of  the 
conducting  boundary: 

-  S7<p  ■  k  =  -Vxy  •  ((c/  +  S)Vxyp)  [12] 

This  is  coupled  with  an  exponential  decay  of  the 
wall  tangential  velocity  with  respect  to  the  core 
flow.  Using  the  divergence  free  condition,  it  is 
possible  to  deduce  all  velocity  components  thus,  if 


desired.  Without  this  integral  treatment,  the 
numerics  were  found  to  be  unstable  for  large 
Hartmann  numbers.  Good  comparisons  were 
obtained  with  respect  to  analytical  solutions. 


(d)  Thick  walls:  This  would  be  the  most  general 
situation,  and  requires  the  solution  to  the  electric 
potential  or  the  magnetic  field  inside  of  the  solid 
wall  region.  While  the  outer  boundary  of  the  walls 
may  be  assumed  to  be  insulating  and  the  conditions 
in  (a),  (b)  be  applied  there,  the  interface  between 
the  walls  and  the  fluid  require  the  continuity  of 
normal  current.  For  an  electric  potential,  this  may 
be  written  simply  in  terms  of  normal  derivatives, 
as: 


dtp 

dn 


=  c  - 


fluid 


dcp 

Jn 


wall 


[13] 


The  electric  potential  may  then  be  solved  inside  the 
wall  from  these  conditions.  In  the  case  of  induced 
magnetic  fields,  the  normal  component  of  the  curl 
of  the  magnetic  field  vector  is  held  constant  across 
the  wall. 


Code  Architecture  and  Attributes 
The  proposed  development  of  a  CFD/MHD  coupled 
solver  will  be  based  on  the  code  architecture  of  an 
existing  unstructured  grid-based  CEM  solver  called 
UPRCS.  This  finite-element  CEM  algorithm  has  been 
implemented  for  massively  parallel  and  scalar  computer 
architectures.  Along  with  the  CEM  solver,  a  complete 
suite  of  user  interface  utilities  are  provided  in  the 
computational  environment  to  go  from  a  CAD  model  to 
the  final  solution.  These  tools  are  shown  in  Figure  [8], 
arranged  as  they  are  used  in  the  simulation  cycle. 

The  tools  include: 

•  The  Rockwell  Science  Center  developed  grid 
generator  UNISG  for  performing  unstructured  grid 
generation. 

•  A  3-D  GUI  (graphical  user  interface)  for  pre  and 
post-processing  called  UNSPREP.  This  tool  aids 
the  user  in  setting  up  grid-related  solver  parameters. 

•  A  help  sensitive  run  control  parameter  input  page 
GUI. 

•  Domain  decompositioning  tools  for  parallel 
simulations. 

•  The  solver,  which  has  been  ported  to  many  parallel 
and  scalar  computer  architectures. 

•  A  collection  of  graphically  driven  post-processing 
tools  for  performing  results  visualization  and  solver 
performance  analysis. 

The  software  written  by  HyPerComp,  is  implemented  in 
a  combination  of  C++,  C,  and  Fortran  and  includes  C++ 
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class  libraries  for  common  re-usable  functions.  It  should 
be  noted  that  all  of  the  software  in  the  solver  suite  is 
either  freely  available  or  developed  in-house. 

Numerical  Schemes:  A  variety  of  numerical  schemes 
have  been  implemented  first  on  structured  meshes  and 
are  being  ported  to  the  unstructured  mesh  framework. 
Schemes  tested  so  far  include  the  Kinetic  Flux  Vector 
Splitting  method  of  Xu  [20],  Roe’s  scheme,  as  presented 
by  Linde  [9],  Runge-Kutta  based  schemes  with  artificial 
dissipation  as  presented  by  Jameson  et  al  [5],  with 
success.  Cases  involving  finite  rate  processes  were 
made  implicit  with  respect  to  the  stiff  source  terms  and 
used  with  a  Roe  scheme.  Details  are  presented  in  [15]. 

Validation 

Due  to  the  inherently  un-intuitive  nature  of  general 
MHD  flows,  validation  poses  a  formidable  problem  in 
code  development  activities.  Some  of  the  earlier  formal 
numerical  work  such  as  that  of  Brio  and  Wu  [2] 
performed  a  rigorous  validation  based  upon  Riemann 
invariants  and  Rankine-Hugoniot  relations,  and  has 
since  been  used  for  validating  other  numerical  schemes. 
Munipalli  et  al.  [15]  have  validated  nonequilibrium  gas 
models  with  MHD  using  analytical  relations  and  power 
balance  estimates  from  Kerrebrock  [6]. 

An  acute  need  exists  for  developing  a  systematic 
validation  procedure  for  MHD  flows  with  various  levels 
of  physics.  A  procedure  based  on  the  method  of 
characteristics  is  being  developed  and  will  be  presented 
in  a  future  paper.  Electrode  surface  phenomena,  electron 
collisional  processes,  ionization  and  recombination  rates 
for  general  methods  of  flow  ionization  all  need  to  be 
well  established  before  introducing  into  modem  code 
environments,  and  this  is  seen  as  a  task  for  the  CFD 
community  at  large  for  the  future.  A  bulk  of  technical 
literature  from  the  1950s  and  1960s  containing 
analytical  and  experimental  results  remains  to  be 
brought  to  bear  in  the  current  context.  A  couple  of 
recent  papers  have  shows  a  firm  and  positive  interest  in 
this  direction  [4],  [17]. 

Sample  Results 

Treatment  of  electric  field  and  current:  As  described 
earlier,  the  electric  field  may  be  computed  from  a 
Poisson-type  equation,  assumed  to  be  zero,  or  backed 
out  from  the  magnetic  field.  Fig.  [2]  shows  a  viscous 
channel  flow  with  a  spatial  distribution  of  magnetic 
field  normal  to  the  plane  of  the  paper,  varying  from  0  to 
1  and  then  back  to  zero.  Current  lines  may  be  shown  to 
lie  in  the  plane  of  the  paper  in  this  case,  and  therefore, 
must  form  closed  loops  or  extend  to  infinity,  in  order  for 
conservation  to  hold.  Contours  of  electric  potential  and 


current  lines  are  shown  here.  Poisson  solutions  are 
obtained  at  each  computational  step. 

Induced  magnetic  field  vectors  are  shown  in  Fig.  [3]  for 
a  supersonic  flow  past  a  ramp,  where  the  wall  normal 
field  (on  an  insulated  wall)  was  assumed  continuous, 
thereby  making  the  induced  field  purely  tangential  to 
the  wall. 

Contours  of  current  density  (computed  from  the  curl  of 
the  magnetic  field)  are  plotted  in  Fig.  [4]  for  a  Mach  4 
flow  past  a  ramp.  The  magnetic  field  is  inclined  along 
the  x-axis,  and  current  lines  extend  to  infinity  in  the  z- 
direction.  Jumps  in  the  current  density  occur  across 
shocks  and  expansion  waves.  A  systematic  validation 
procedure  by  which  these  jumps  may  be  verified  is 
being  developed. 

Inlet  flow  fields:  Flow  fields  in  various  inlets  designed 
to  operate  at  Mach  numbers  between  6  to  10  have  been 
studied.  Sample  results  in  Figure  [5]  show  a  strong 
dependence  of  the  induced  magnetic  field  on  the  flow 
structure.  An  induced  magnetic  field  about  10  times 
larger  than  the  applied  field  in  the  x-direction  has  been 
observed  for  a  perfectly  conducting  gas.  When  finite 
values  of  electrical  conductivity  is  used,  the  induced 
field  is  diminished  greatly,  and  a  more  complete 
correlation  using  the  Rankine-Hugoniot  relations  is 
being  currently  derived. 

Figure  [6]  shows  the  increase  in  the  shock  angle  with 
increased  applied  magnetic  field  in  the  x-direction.  A 
small  discontinuity  developing  close  to  the  wall  is  the 
second  magnetosonic  shock  whose  angle  increases  with 
applied  field.  A  peak  pressure  value  near  the  wall 
corresponding  to  this  second  shock  is  seen  to  increase 
with  applied  field.  Figure  [7]  shows  an  implementation 
of  MHD  as  a  source  term  in  an  Euler  solver  for  the  flow 
past  a  ramp.  Pressure  profile  along  the  exit  plane  shows 
a  linear  increase  in  pressure  (corresponding  to  the 
pressure  gradient  term,)  balanced  by  the  Lorentz  force. 
The  numerical  value  of  this  gradient  may  readily  be 
validated. 

An  unstructured  grid  implementation  of  these  concepts 
has  been  initiated.  Figure  [1 1]  shows  a  3-D  unstructured 
flow  in  a  Mach  8  inlet  with  and  without  an  applied  field. 
A  weak  Mach  wave  along  the  side  of  the  inlet  was  seen 
to  intensify  to  a  shock  upon  application  of  a  magnetic 
field.  The  shock  angle  at  the  inlet  is  also  altered,  and  a 
greater  compression  is  achieved  for  the  same  flow 
turning  angle. 

Viscous  effects:  Hypersonic  flows  on  slender  fore¬ 
bodies  often  involve  the  development  of  thick  boundary 
layers,  which  interact  with  the  shock  waves  to  form 
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extremely  complex  shock  layer  dynamics.  Sample 
results  from  these  with  and  without  MHD  effects  are 
shown  in  Figs.  [9],[10].  A  more  detailed  study  of 
viscous  effects  with  induced  fields  is  underway. 
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Magnetic  Field  Wall 


Figure  1  :  Electric  field  in  a  channel  with  segmented  electrodes 


Figure  2  :  Electric  current 
vectors  (above)  and  electric 
potential  contours(below)  for 
MHD  flow  in  a  channel  where 
an  out  of  plane  magnetic  field 
is  applied,  ramping  up  between 
x=l,1.5  and  down  between 
2.5,3. 
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Induced  Magnetic  Field 


Figure  3  :  Induced  magnetic  field  vectors  for 
ramp  flow  problem  with  an  imposed  vertical 
field. 
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Figure  4:  Contours  of  electric 
current  derived  from  curl  (B)  for 
Mach  4  flow  past  a  ramp.  Current 
travels  in  the  z-direction  to  oo. 
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Figure  5:  Mach  8  inlet.  Clockwise  from  top-left:  (a)  Density  contours  with  B  =  0;  (b),(c)  x,y 
components  of  the  magnetic  field;  (d)  Magnitude  of  the  induced  field  (-0.3  to  10  *  Bx  applied) 


8 

American  Inslimte  of  Aeronautics  and  Astronautics 


P(Pa) 


DCC?  1  AVAIL  ABLE  copy 


Figure  6:  Increasing  magnetic  field  strength  (0.02  to  0.08  Tesla)  on  a  Mach  4  flow  (ideal  MHD)  solution. 
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Figure  7:  Pressure  along  the  exit  plane  of  a  Mach  4  flow  past  a  ramp  with 
MHD  added  as  source  terms.  Perfect  gas  post-shock  pressure  has  been  marked 
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Figure  8:  UPRCS  solver  code  suite  developed  by  HyPerComp,  Inc. 


Figure  9:  Mach  10  inlet  (from  inviscid  design)  with  viscous  effects  and  large  mass  spillage 
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Figure  10  :  Inlet  from  Fig.  [9]  with  applied  magnetic  field 
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Figure  1 1  :  Preliminary  unstructured  mesh  simulations  of  hypersonic  inlet  with  MHD  interactions. 
Above  left:  B  =  0,  below  left:  Bx  -  1  Tesla,  a  =  40  IQIm  in  a  Mach  8  inlet.  Electrical  conductivity  a 
siirmle  temperature  dependant  function. 
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Abstract 

A  numerical  platform  based  on  detailed  physics 
for  the  simulation  of  plasma/MHD  aerodynamic 
flows  is  under  development  at  HyPerComp,  Inc. 
In  our  earlier  papers  on  the  subject  ([l]-[4]),  we 
have  addressed  the  treatment  of  seeded  and 
unseeded  air  using  a  variety  of  models  to 
describe  the  electromagnetic  field.  These  studies 
have  been  migrated  to  an  unstructured  mesh 
based  parallel  code  structure,  which  is  being 
currently  used  to  study  MHD  effects  on 
hypersonic  inlets. 

We  present  in  this  paper,  our  ongoing  studies  of 
ionization  mechanisms  induced  by  e-beams. 
The  principal  difference  here  will  be  the 
formation  of  ions  at  low  temperature  caused  by 
the  electron  beams,  as  opposed  to  the  high 
temperature  formation  of  NO  and  its  ions  (or  the 
equilibrium  ionization  of  alkali  salts)  in  our 
previous  studies.  Computations  were  performed 
on  a  typical  3-D  supersonic  inlet  geometry.  E- 
beam  induced  plasma  kinetics  were  coupled 
with  MHD  equations  and  3D  inviscid  flow 
equations.  A  two-temperature  (T,  Tv/Te)  gas 
model  is  used  to  include  effects  of  thermal  non¬ 
equilibrium.  Effects  of  e-beam  induced 
ionization  and  MHD  forces  on  supersonic  flow 
deceleration  were  investigated. 

Introduction 

Weakly  ionized  gases  have  a  variety  of 
aerospace  applications.  MHD  inlet  control, 
supersonic  drag  reduction,  plasma-enhanced 
combustion  (in  scramjets)  are  some  of  the 
notable  examples.  Simulation  of  such  flows 
requires  the  solution  of  Navier-Stokes  equations, 
finite-rate  chemical  kinetics  and  Maxwell’s 
equations.  Frequently,  these  flows  occur  in 
complex  3D  geometries  and  hence  can  be  very 
computationally  intensive.  Higher  order 


numerical  schemes  are  required  in  order  to 
accurately  capture  the  shock  formation  and 
reflection  phenomena.  Hence  a  design  tool  to 
simulate  such  flows  would  have  to  use  high- 
performance  computational  methods  with  robust 
higher  order  numerical  schemes  to  obtain 
solutions  in  a  reasonable  amount  of  time. 

High  Speed  Inlet  control  using  MHD  and 
plasma-assisted  combustion  of  complex 
hydrocarbon  fuels  in  scramjets  are  applications 
where  such  a  design  tool  would  be  invaluable.  A 
key  issue  in  these  applications  is  ionizing  the 
incoming  air  stream.  The  requirements  for 
producing  this  ionized  air  stream  are  very 
stringent.  Large  volumes  (~  1  m3),  low  gas 
temperatures  (1000  -  2000  K),  high  electron 
density  (~  1013  /cm3),  long  residence  times  (>  10 
msec),  low  power  budget  (1  -  10  MW/m3), 
strong  magnetic  fields  (~  10  Tesla)  are  the  main 
requirements  for  these  air  plasmas.  Several 
techniques  have  been  proposed  in  recent 
literature,  to  achieve  these  conditions.  Some  of 
these  techniques  are,  high-energy  electron 
beams  (10  <  E  <  100  KeV),  high  voltage 
nanosecond  pulses,  DC  &  RF  discharges,  E- 
beams  in  vibrationally  excited  gases  and  seeding 
of  the  incoming  gas  stream.  High-energy 
electron  beams  are  most  efficient  in  terms  of 
ionization  efficiency  [5]  and  hence  have  been 
used  in  several  studies  to  model  air  plasmas. 

Electron  beams  produce  directed  ionization  in 
gases.  They  provide  a  convenient  means  to 
generate  locally  high  electrical  conductivity  at 
relatively  low  temperatures.  The  energy  of  the 
e-beams  can  be  varied  in  the  range  of  10  to  100 
keV,  to  provide  high  enough  electrical 
conductivity  at  moderate  flow  static  pressures. 
The  efficiency  of  electron  beams  diminishes 
with  increasing  pressure  due  to  the  inherent 
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nature  of  the  electrons  to  be  lost  with  greater 
ease  in  such  situations  (largely  due  to  attachment 
with  02).  The  e-beam  chemistry  model  used  in 
this  work  has  been  used  by  several  authors  [6-8] 
and  has  been  validated  by  experiment  [8].  The 
reactions  and  rate  constants  are  those  used  in  [8] 
and  are  shown  in  Table  1  in  the  Appendix. 

The  tools  being  developed  here  fill  the  need  for 
a  high-performance,  numerically  robust  design 
capability  to  accurately  predict  flow  fields  of 
weakly  ionized  gases.  Such  a  design  tool  would 
include  real  gas  effects  such  as  variable  transport 
properties,  finite  electrical  conductivity,  finite- 
rate  chemical  kinetics,  thermal  non-equilibrium 
and  effects  of  electromagnetic  fields. 

Numerical  formulations  and  perspectives 


Maxwell  -  Navier-Stokes  coupled  equation  set 


The  conservation  equations  governing  the  flow 
of  weakly  ionized  gases  may  be  written  in  the 
following  form: 

—  +  V  •  (Fx,  Fy,  Fz )  +V  •  (Hx,  Hy,  Hz)  =  S 
dt 

in  which  the  fluxes  F  are  inviscid  or  convective, 
while  H  are  viscous  or  dissipative.  The  source 
term  S  comprises  of  the  production  of,  and 
interaction  between  the  various  flow  quantities. 
This  equation  may  be  expressed  in  the  integral 
conservation  form  in  a  control  volume  Q 
bounded  by  8Q,  as  follows: 


J (j  [F„ds  +  H„ds]=  JSc/Q 

n  so.  n 

Fn  and  H„  in  the  above  equation  are  face  normal 
flux  components.  The  vector  used  in  this 
equation  are  expanded  below.  An  index  s  is  used 
for  the  various  gaseous  species  that  are 
considered.  The  gas  model  used  in  this  work 
comprises  of  the  following  components:  N2,  02, 
N2+,  02+,  02,  e.  We  currently  employ  a  two 
temperature  model  of  ionizing  air,  similar  to  that 
of  Park  [9].  Energy  addition  by  MHD  at  high 
temperatures  can  bring  about  a  departure  from 
equilibrium  between  the  various  energy  modes. 
There  are  two  energy  equations,  one  for  the  total 
energy  conservation  and  the  other  describing  the 


vibrational-electronic  energy.  The  following 
matrices  describe  the  conservation  laws: 


V  =  [u,v,w] 
PU 

r  s  n 

_  pUVT+pn 
PUH 
pU  e 

n  v 


-  p°?ys  ■ « 

-  r 

n 

V  ■  f  -  rjVT  ■  n  -  tjVT  ■  n  -  hP?ys 

-TjVT-n-p^P'PPy," 


0  ><  2) 

J  •  E'  +  6)  +2 


e,rad 


The  transfer  of  energy  across  the  various  energy 
levels  is  contained  in  the  nonequilibrium  energy 
production  rate^yv.  This  term  contains  the 
energy  transfer  to  the  electronic  and  vibrational 
modes  via  collisions,  electron  pressure  gradient, 
Landau-Teller  type  of  relaxation  of  these  modes 
with  the  translational  and  rotational  modes,  and 
release  in  vibrational  and  electronic  energy 
caused  by  ionization  or  the  change  in  the 
number  densities  of  polyatomic  molecules.  The 
expressions  used  to  model  these  features  may  be 
found  in  various  references  (e.g.:  Munipalli  [1- 

4]). 

Radiative  energy  addition  by  means  of 
microwave  and  other  beams  can  possibly  be 
used  to  cause  an  increase  in  the  electron  energy 
and  ionization.  This  is  modeled  by  the  terms 
Qrad.  Fluxes  Fn  and  Hn  are  written  at  the  cell 
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face  in  each  control  volume.  In  the  results 
presented  in  this  paper  effects  of  dissipative 
fluxes  (i.e.  viscosity,  thermal  conductivity, 
species  diffusivity  and  mobility  effects)  are  not 
included. 

MHD  Model: 


In  the  calculations  presented  in  this  paper,  the 
magnetic  Reynolds  number  is  assumed  to  be 
small.  This  implies  that  the  induced  magnetic 
field  is  negligible  and  that  the  electric  field  may 
be  derived  from  a  scalar  potential:  E-  -V<p  ■ 
Ohm’s  law  may  then  be  used  to  relate  this 
potential  to  the  current  density: 

J  =  <t(-V<p  +  VxB)-— (jxB) 

B 

Where  a  denotes  the  electrical  conductivity, 
(ox/B  is  the  Hall  parameter,  J  is  the  electric 
current  density  and  B  is  the  applied  magnetic 
field  vector.  Since  the  induced  magnetic  field  is 
neglected,  the  only  electromagnetic  field 
quantity  that  needs  to  be  computed  numerically 
at  each  time  step  is  the  electric  potential. 

When  the  divergence  of  the  above  expression  is 
set  to  zero  and  is  re-written  in  terms  of  the 
electric  potential,  a  Poisson  type  equation  is 
obtained,  which  must  be  solved  at  each 
computational  time  step.  In  the  results  presented 
in  this  paper,  the  Hall-effect  terms  are  neglected 
in  the  generalized  Ohm’s  law  given  above. 
Setting  the  divergence  of  the  current  density 
(without  the  Hall  terms)  to  zero  leads  to  the 
following  expression: 

V-[<r(-Vp  +  Kx5)]=0 
Solution  of  the  above  Poisson  equation  yields 
the  electric  potential  (p  and  electric  current 
density  j.  The  joule  heating  term  and  Lorentz 
forces  in  the  energy  and  momentum  equations 
can  thus  be  obtained.  In  this  work,  Bx  =  Bz  =  1 
T,  By  =  0. 


Chemistry  model: 

The  chemistry  model  used  in  this  work  is  based 
on  Ref  [8].  This  model  has  been  validated  with 
experimental  data  and  hence  has  been  chosen  in 
this  work.  The  electron  production  rate  by  the  e- 
beam  (S0)  is  chosen  to  be  constant  and  set  equal 


to  1.5E24  m'3s'\  The  spatial  extent  of  the  e- 
beam  is  set  by  assuming  (see  Figure  1) 

S  =  S0,  for  X0  <  X  <  Xi 


S  =  0,  for  X  <  Xo  or  X  >  Xi 


The  electrical  conductivity  of  the  ionized  air  is 
computed  using  the  expression  for  partially 
ionized  gas  given  by  Rosa  [10] 


<j  = 


nee 


mece 


^"kQk  +  3-9/J; 


<  e2  ^ 
8  7ce0kT  j 


log  A 


where  A  = 


1.24*10  i 


,7  T1.5  \ 


n 


0.5 


In  a  mixture  of  gases,  n,  denotes  the  number 
density  of  ions  and  Qk  is  the  neutral-neutral 
collision  cross-section. 


Boundary  Conditions: 

Boundary  and  initial  conditions  for  the 
governing  equations  described  above  must  be 
specified  accurately  in  order  to  prevent 
numerical  instabilities  and  non-physical 
behavior  of  the  solution.  For  supersonic  and 
super-alfVenic  flows,  information  propagates 
only  in  the  downstream  direction.  The  inflow 
variables  are  then  fixed  and  outflow  variables 
are  obtained  by  extrapolation.  For  subsonic  or 
sub-alfvenic  flows,  characteristic  relations  must 
be  used  in  order  to  propagate  information  either 
into  or  out  of  the  domain. 

No  slip  boundary  conditions  are  applied  to  the 
velocity  along  the  wall  for  viscous  flows.  In  the 
present  work,  the  flow  is  considered  to  be 
inviscid,  hence  the  gradient  of  the  tangential 
component  of  velocity  is  set  equal  to  zero  at  the 
walls.  Normal  derivatives  for  species  densities 
are  set  equal  to  zero  on  the  walls.  The  gradient 
of  the  gas  temperature  and  vibrational 
temperature  is  set  equal  to  zero  at  the  walls. 
Electromagnetic  boundary  conditions  may  be 
more  involved,  depending  on  the  choice  of  wall 
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material.  In  the  present  work,  the  electric 
current  density  normal  to  all  boundaries  is  set 
equal  to  zero,  i.e. 


Jn  =J.n  =  0 


Snace/time  integration:  Numerical  Scheme 


A  point-implicit  finite  volume  scheme  is  used  to 
solve  the  governing  equations  shown  earlier.  For 
a  cell  bearing  the  index  i,  volume  Q„  the  time 
advancement  scheme  is  written  as: 


U'+A'  =U‘  +C~xr‘ 

where, 


CL  = 

and 


=  M\m  +  M\isc  -  Q,M'src 

The  right  hand  side  vector  r  is  the  flux 
summation  based  on  quantities  at  time  level  t. 
Inviscid  fluxes  are  computed  from  a  Roe-type 
procedure,  in  which  a  density  based  averaging  is 
performed  similar  to  the  perfect  gas  equivalent. 
Viscous  and  dissipative  fluxes  are  computed 
from  central  differences.  Details  of  this 
implementation  may  be  found  in  Refs.  [4,12]. 


In  the  above  equation,  /  is  an  n  x  n  (n  is  the  total 
number  of  dependent  variables  in  the  vector  U) 
identity  matrix  and  ML’  is  the  combined 
Jacobian  including  inviscid,  viscous  and  source 
term  contributions. 


Reacting  and  plasma  flow  equations  are 
characterized  by  widely  varying  spatial  and 
time-scales.  The  point-implicit  formulation  is 
necessary  in  order  to  obtain  stable  solutions  for  a 
stiff  system  of  governing  equations.  As 
mentioned  earlier,  Mvisc  =  0  in  the  results 
presented  in  this  paper. 

The  formal  order  of  accuracy  of  this  numerical 
procedure  is  enhanced  by  using  the  TVD 
approach,  as  described  by  Barth  [11],  using  local 
gradients  computed  from  a  Gaussian  summation. 
Work  is  underway  to  develop  this  environment 


into  a  higher  order  multi-physics  solver  using 
the  Discontinuous  Galerkin  technique. 

Results 

Figure  1  shows  the  geometry  used  in  the  present 
work.  The  e-beam  chemistry  model  was 
validated  by  conducting  a  simulation  with  free 
stream  conditions  of  1  atm,  300  K  (as  in  Ref 
[8]).  For  the  purpose  of  validation,  the  electron- 
ion  formation  by  means  of  the  e-beam  is 
simulated  in  the  region  3.46  <  x  <  4.56m.  The 
production  rate  for  electron-ion  pairs  by  the  e- 
beam  is  assumed  to  be  1.5E24  m'3/s  as  in  [8]. 
The  magnitude  of  the  electron  concentration  is 
compared  with  results  reported  in  Fig.  9  in  ref. 
[8].  Figure  2  shows  the  axial  variation  of 
electron  number  density  for  the  above 
conditions.  It  is  see  that  the  peak  electron 
number  density  is  close  to  that  predicted  by 
experiments  (1.8E16  #/m3-s)  in  Ref  8. 

We  use  a  base  inlet  geometry  designed  for  Mach 
8  operation  from  inviscid  considerations,  earlier 
studied  by  Sheikin  et  al  [13].  A  flight  Mach 
number  of  6  at  an  altitude  of  30  km  was  chosen 
as  an  off-design  condition  in  our  study. 
Corresponding  to  this  altitude,  the  ambient 
density,  pressure  and  temperature  are  1.84E-2 
Kg/m3,  1.2E3  N/m2  and  226.5  K,  respectively. 
With  an  intent  to  improving  the  performance  of 
this  inlet  at  the  lower  than  design  Mach  number, 
the  following  case-studies  were  made. 

Flow  Without  MHD 

Figure  3(a),  (b)  show  contours  of  axial  velocity, 
and  temperature.  The  velocity  and  temperature 
contours  show  the  formation  of  shocks  due  to 
the  two  ramps.  Shock  reflection  from  the  top 
boundary  (inviscid  wall)  is  also  seen.  Figure  4 
show  the  axial  variation  of  gas  temperature  and 
vibrational  temperature  along  the  ramp  surface. 
It  is  seen  that  there  is  very  little  change  in  the 
vibrational  temperature  of  the  gas  along  the 
length  of  the  ramp  channel.  At  the 
pressure/temperature  conditions  being 

considered  here,  the  V-T  relaxation  time  is 
rather  large  (a  few  milli-seconds,)  and  Tv  lags 
far  behind  the  translational  temperature  T.  When 
Joule  heating  and  other  forms  of  energy  addition 
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preferentially  excite  the  vibrational/electronic 
modes,  the  vibrational  temperature  can  grow  to 
be  in  excess  of  the  translational  temperature,  and 
energy  is  exchanged  across  these  modes  by 
relaxation.  This  may  be  seen  in  the  subsequent 
results  in  this  section. 

Constant  conductivity  with  Bx  =  Bz  =  1: 

A  parametric  study  of  the  effect  of  electrical 
conductivity  and  extent  of  the  ionized  region 
was  also  performed.  These  simulations  were 
made  on  a  coarse  mesh.  Figure  5  shows  the 
effect  of  the  MHD  forces  on  the  flow  field.  For 
this  case  the  electrical  conductivity  is 
maintained  at  20  mho/m  in  a  region  between  1 .0 
<  x  <  5.5  m.  Figure  5(a)  and  (b)  show  the 
velocity  contours  with  and  without  MHD 
effects.  It  is  seen  that  there  is  a  noticeable 
change  in  shock  angle  and  flow  deceleration  due 
to  the  MHD  forces. 

Constant  conductivity  with  Bx  =  0,  Bz  =  1: 

Figure  6  shows  the  current  lines  for  a  case  with  a 
constant  conductivity  of  20  mho/m  but  with  a 
constant  magnetic  field  of  Bz  =  1 .  The  absence 
of  an  X-component  of  the  magnetic  field  and  the 
insulating  boundary  conditions  give  rise  to 
current  lines  that  close  on  themselves  in  the  XY 
plane. 

MHD  with  variable  a  and  e-beam  ionization 

MHD  forces  generated  in  ionized  air  depend 
primarily  on  the  electrical  conductivity  and  the 
applied  magnetic  field.  The  electrical 
conductivity  of  air  is  governed  by  the  number 
density  of  mobile  charge  carriers  in  the  flow  (see 
expression  for  electrical  conductivity  described 
earlier).  Hence,  it  is  instructive  to  estimate  the 
conductivity  of  e-beam  ionized  air  as  a  function 
of  S0  (e-beam  induced  ionization). 

A  few  parametric  cases  were  run  to  study  the 
effect  of  beam-induced  electron  production  rate 
(S0)  on  the  maximum  electrical  conductivity. 
These  values  are  reported  in  Table  2  in  the 
appendix.  It  is  known  that  the  peak  electron 
density  is  proportional  to  the  square  root  of  the 
e-beam  induced  electron  production  rate  [7])  and 


this  is  seen  in  Table  2.  The  electron  production 
rate  S0  is  determined  by  the  e-beam 
characteristics  and  is  given  by  the  following 
expression  [7]: 


where,  W;  is  the  ionization  energy  to  produce  an 
electron-ion  pair  (34  eV),  Eb  is  the  energy  of  the 
e-beam  (typically  ~  100  KeV),  p  is  the  gas 
density  (and  Y  is  the  electron  stopping  power  in 
(MeV  m2/kg).  At  higher  altitudes,  the  gas 
density  is  lower,  hence  for  a  given  e-beam 
energy,  a  proportionally  higher  beam  current 
density  is  required  to  produce  the  required  e- 
beam  induced  electron-ion  production  rate.  The 
electron  production  rate  of  E24  m'3/s  in  Ref  [8] 
is  based  on  a  beam  current  density  of 
approximately  5  mA/cm2  and  beam  energy  of 
about  30  KeV.  Column  II  in  Table  2  shows  the 
current  density  scaled  with  gas  density  for  the 
various  beam  induced  production  rates  shown  in 
column  I.  It  is  seen  that  to  produce  an  electrical 
conductivity  around  20  mho/m,  a  very  high 
value  of  beam  current  is  required  and  hence 
could  present  practical  engineering  problems  in 
implementation.  For  the  beam  strength  of  Sc  = 
1.5E24  m'3/s  used  in  this  work,  the  maximum 
value  of  sigma  is  about  5  mho/m. 


Figure  7  (a),  (b),  (c)  show  contours  of  axial 
velocity,  temperature  and  electron  number 
density,  respectively.  The  free-stream  Mach 
number  is  6,  and  the  e-beam  is  applied  in  the 
region  from  4.0  <  x  <  4.75  m.  The  strength  of 
the  e-beam  induced  ionization  (S0)  was  equal  to 
1.5E24  m'3/s.  A  uniform  magnetic  field  of  Bx  = 
Bz  =  IT  was  used  in  these  simulations.  The 
velocity  and  temperature  profiles  look  similar  to 
the  case  where  there  was  no  MHD  present.  It  is 
seen  in  Fig  7(c)  that  in  the  region  with  the  e- 
beam  the  electron  density  is  high  but  close  to 
zero  everywhere  else  in  the  domain.  Figure  8 
(a)  shows  the  current  vectors  in  a  cross-section 
at  x  =  4.37  (center  of  the  e-beam  region). 
Figure  8  (b)  shows  the  current  lines  with 
contours  of  cr.  The  current  traces  are  spiraling 
due  to  the  presence  of  a  magnetic  field  in  the  X- 
direction. 
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Figures  9  (a)  and  9  (b)  show  the  variation  of 
axial  velocity  and  vibrational  temperature  near 
the  ramp  surface,  respectively.  Figure  9  (a) 
show  the  reduction  in  the  axial  velocity  due  to 
the  presence  of  the  JXB  terms.  It  is  seen  in  Fig. 
9  (b)  that  addition  of  the  Ohmic  heating  terms 
(J.E)  quickly  increases  the  vibrational 
temperature  as  compared  to  the  case  with  no 
MHD. 

Variable  a,  S0  =  1.5E26  ni"3/s  and  e-beam 
applied  over  3.5  <  x  <  5.5 

Figure  10  (a)  shows  variation  of  T,  Tv  and 
electron  number  density  with  S0  =  1 .5E26  m'3/s. 
Figure  10  (b)  shows  variation  of  T,  Tv  and 
electrical  conductivity.  The  higher  value  of 
beam  induced  ionization  gives  rise  to  a  higher 
value  of  electrical  conductivity  and  electron 
number  density.  Since  the  contribution  of  the 
Ohmic  heating  term  to  the  vibrational  mode  is 
high,  the  vibrational  temperature  is  high  as 
pointed  out  above. 

Conclusion 

This  paper  presents  the  results  of  our  ongoing 
investigations  in  the  inclusion  of  real  gas  effects 
in  a  practical  flow  solver  to  simulate  plasma  and 
MHD  effects  in  hypersonic  flow.  The  cold 
incoming  supersonic  stream  is  ionized  by  the  use 
of  e-beams.  A  plasma  kinetic  model  describing 
the  ionization  of  cold  air  is  coupled  with  a 
Poisson  solver  and  3-D  MHD  equations.  It  is 
seen  that  large  volume  E-beam  generated 
ionization  is  required  to  bring  about  reasonable 
amounts  of  shock  displacement.  Inviscid  flow 
simulations  show  that  high  beam  current 
densities  are  required  to  generate  the  ionization 
levels  needed  for  shock  displacement.  Real  gas 
effects  such  as  species  diffusion,  mobility  and 
viscosity  will  be  included  in  future  simulations. 
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Appendix 


Table  1 :  E-beam  chemistry  model  with  reaction  rates  constants 


No. 

Reactions 

A 

n 

Ea/R 

1 

e'  +  202  =  02'+  02 

2.5E-42  (m6/s) 

0 

0 

2 

e"  +  O2+  N2  =  O2"  +  N2 

1.6E-43  (m6/s) 

0 

0 

3 

e  +  O2  =  O2 

2.0E-12  (m3/s) 

0 

0 

4 

e'  +  N2+  =  N2 

2.0E-12  (m3/s) 

0 

0 

5 

O2  +  O2  =  2O2  +  e’ 

2.2E-24  (m3/s) 

0 

0 

6 

O2”  +  N2  =  O2  +  N2  +  e’ 

1.8e-26  (m3/s) 

0 

0 

7 

O2  +  02  +  O2  =  3  02 

1.55E-37  (m6/s) 

0 

0 

8 

O2’  +  02+  +  N2  =  2  O2  +  N2 

1.55E-37  (m6/s) 

0 

0 

9 

O2’  +  N2+  +  O2 =  2  O2  +  N2 

1.55E-37  (m6/s) 

0 

0 

10 

O2"  +  02+  +  N2  —  O2  +  2N2 

1.55E-37  (m6/s) 

0 

0 

Table  2:  Variation  of  Electrical  conductivity  with  pressure  and  S0 


So  (#/m3/s) 

Jb  (m  A/cm2) 

p(kg/m3) 

Pressure 

(N/m2) 

Ne(#/m3) 

a  (mho/m) 

1.5E24 

5 

1.22 

1.01E5 

1.8E16 

0.0002 

1.5E24 

3.32E2 

1.84E-2 

1.2E3 

7.9E17 

5.5 

1.5E25 

3.32E3 

1.84E-2 

1.2E3 

2.64E18 

11.7 

1.5E26 

3.32E4 

1.84E-2 

1.2E3 

8.53E18 

19.56 

1.5E23 

3.32E1 

1.2E3 

2.6E17 

2.5 

1.5E22 

3.32 

1.84E-2 

1.2E3 

7.7E16 

0.8 
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Figure  1 :  Schematics  of  the  geometry  (all  dimensions  in  meters) 
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Fig  7:  Axial  velocity  (a),  temperature  (b)  and  electron  number  density  (c)  with  MHD 


(a)  (b) 

Figure  8:  Current  vectors  (a)  and  current  lines  (b)  in  the  e-beam  region 
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(a)  (b) 

Figure  9:  Axial  variation  of  velocity  (a)  and  vibrational  temperature  (b)  with  MHD 


Figure  10:  Axial  variation  of  T,  Tv  and  electron  number  density  (a),  and  T,  Tv  and  a  (b)  for  a  case  with  S0 

=  1 .5E26  m'3/s 
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Effect  of  wall  conduction  in  proposed  MHD  enhanced  hypersonic  vehicles 

Ramakanth  Munipalli+,  Shashi  Aithal’  and  Vijaya  Shankar* 
HyPerComp,  Inc.,  31255,  Cedar  Valley  Dr.,  Suite  #327, 

Westlake  Village,  CA  91362 


Abstract 

In  AJAX-type  vehicle  design  concepts  (e.g.,  ref. 
[7]),  the  flow  of  a  conducting  medium  in  the 
presence  of  a  magnetic  field  produces  currents 
that  must  form  closed  loops  or  extend 
indefinitely,  in  the  absence  of  charge  sources. 
While  the  solid  walls  of  the  airframe  are 
normally  assumed  to  be  insulating  in  CFD 
estimates  of  AJAX  performance,  even  mildly 
conducting  airframe  materials  cause  currents  to 
leak  away  from  the  flow  (considering  the  very 
low  achievable  values  of  electrical  conductivity 
of  air  in  flight,)  and  may  potentially  cause  major 
changes  in  flow  compression  and  shock 
structures  on  this  account.  Here,  we  present  the 
effect  of  conducting  walls  in  hypersonic  inlet 
flow  computations.  Even  small  cracks  in 
insulations  (4-5  orders  of  magnitude  smaller 
than  the  flow  length  scale,)  have  large  scale 
effects  on  the  flow,  and  the  influence  of  MHD 
on  flow  control  aspects.  Conducting  wall 
physics  is  superposed  on  our  previous  work  in 
this  area  (refs.  [2], [13]-[17]). 

Introduction 

Recent  years  have  witnessed  a  great  surge  in 
simulation  capabilities  pertaining  to  ionizing  air 
flows  and  MHD  (e.g.,  ref.  [20]).  This  has  been 
coupled  with  experiments  demonstrating 
physical  phenomena  of  interest  in  such  flows, 
and  analytical  estimates  of  performance 
enhancements  in  hypersonic  vehicle  design 
caused  by  these  techniques. 

Effect  of  wall  conductivity 
Existing  literature  in  the  area  of  hypersonic 
flows  with  MHD  frequently  relies  on  the 
assumption  that  the  walls  bounding  air  flows  are 
insulating,  or  have  some  pre-set  potential 
distribution  that  determines  the  current  within 
the  fluid  medium.  When  the  walls  are  fully 


conducting,  or  comprise  of  insulation  in  which 
imperfections  may  occur  (albeit  minute,) 
literature  has  shown  numerous  situations  in 
which  the  computed  pressure  is  vastly  different 
(by  more  than  an  order  of  magnitude,)  from  the 
insulating  wall  assumption.  (See,  e.g.,  Refs 
[4], [11]).  In  airframes,  a  metallic  structure  is 
more  of  a  norm  than  an  exception,  and  this 
seems  to  call  for  a  major  revision  of  MHD 
performance  estimates  obtained  from  insulating 
wall  approximations. 

Figure  [1]  shows  a  situation  in  which  current  is 
generated  in  a  flow  by  virtue  of  conducting  gas 
moving  in  the  presence  of  a  magnetic  field.  This 
current  does  not  penetrate  insulating  walls 
(shown  in  thick  black,)  but  seeks  out  the  high 
conduction  path  through  cracks  in  insulation  or 
solid  metallic  elements  in  order  to  close  loops. 
In  this  paper,  we  attempt  to  assess  the  MFID 
effects  of  conducting  walls  and  quantify  them 
for  cases  of  interest.  The  pressure  gradient  is 
accentuated  (10s  of  times)  when  more  current 
flows  in  the  fluid,  typically,  when  conducting 
walls  are  present.  Indeed,  any  MHD  effect  that 
is  felt  upon  a  flow  is  dependent  on  the  electrical 
conductivity  of  the  gas.  A  practical  realization  of 
these  concepts  in  a  flight  vehicle  will  depend 
critically  on  enhanced  electrical  conductivity 
and  the  ability  to  carry  high  strength  magnets 
on-board. 

Ionization  and  electrical  conductivity  of  gases 

The  electrical  conductivity  of  a  gaseous  medium 
such  as  air  may  be  enhanced  by  various  forms  of 
ionization.  Weakly  ionized  gases  have  a  variety 
of  aerospace  applications.  MHD  inlet  control, 
supersonic  drag  reduction,  plasma-enhanced 
combustion  (in  scramjets),  are  some  of  the 
notable  examples.  A  key  issue  in  these 
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applications  is  ionizing  the  incoming  air  stream. 
The  requirements  for  producing  this  ionized  air 
stream  are  very  stringent.  Large  volumes  (~  1 
m3),  low  gas  temperatures  (~  2000  K),  high 
electron  density  (~  1013  /cm3),  long  residence 
times  (>  10  msec),  low  power  budget  (1  -  10 
MW/m3),  strong  magnetic  fields  (~  10  Tesla)  are 
the  main  requirements  for  these  air  plasmas. 

Several  techniques  have  been  proposed  in  recent 
literature,  to  achieve  these  conditions.  Some  of 
these  techniques  are:  high-energy  electron 
beams  (E  >  100  KeV),  high  voltage  nanosecond 
pulses,  DC  &  RF  discharges.  E-beams  in 
vibrationally  excited  gases  and  seeding  of  the 
incoming  gas  stream.  We  consider  high  energy 
electron-beams  to  produce  ionization  in  air  in 
this  paper.  Some  recent  studies  have  indicated 
their  superiority  to  other  methods  of  ionization 
(ref.  [9]). 

Electron  beams  produce  directed  ionization  in 
gases.  They  provide  convenient  means  to 
generate  locally  high  electrical  conductivity  at 
relatively  low  flow  temperatures.  Ref[2] 
presents  our  implementation  of  the  electron 
beam  models  in  a  thermochemical 
nonequilibrium  MHD  solver. 


Problem  Formulation 

Flow  Equations 

Over  the  years,  we  have  developed  models  of 
varying  complexity  to  describe  MHD  effects  on 
air.  These  models  range  from  ideal  gas,  to 
equilibrium  air  with  seeding,  nonequilibrium 
ionizing  air  with  multiple  temperatures,  and 
have  recently  extended  these  to  include  the 
effects  of  electron  beams.  The  reader  is  referred 
to  Munipalli  et  al.  [13-17]  for  a  description  of 
the  electromagnetic  and  gas  models  and  our 
previous  experience  with  MHD  modeling  in 
hypersonic  flow  acceleration  and  inlet  studies. 

Electromagnetics  -  the  inductionless  approach 
In  our  earlier  work  (refs.  [1 3]-[l 7]),  we 
considered  the  differences  between  various 
electromagnetics  models  for  use  in  magneto¬ 
aerodynamics,  ranging  from  fully  coupled  NS- 
Maxwell  equations  to  various  forms  of 


inductionless  models  using  the  electric  potential 
(refs,  [14],  [16]).  Electrical  conductivity  of  air, 
with  most  practical  ionization  techniques  being 
investigated  in  the  present  day,  does  not  exceed 
10-100  mho/m,  resulting  in  a  magnetic  Reynolds 
number  that  is  extremely  small.  In  the 
calculations  presented  in  this  paper,  the 
magnetic  Reynolds  number  is  assumed  to  be 
small.  This  implies  that  the  induced  magnetic 
field  is  negligible  and  that  the  electric  field  may 
be  derived  from  a  scalar  potential:  E  =  -V<p. 
Ohm’s  law  may  then  be  used  to  relate  this 
potential  to  the  current  density: 

J  =  <r(-V<p  +  VxS)-—  (jxB) 

B 


Where  a  denotes  the  electrical  conductivity, 
cot/B  is  the  Hall  parameter,  J  is  the  electric 
current  density  and  B  is  the  applied  magnetic 
field  vector.  Since  the  induced  magnetic  field  is 
neglected,  the  only  electromagnetic  field 
quantity  that  needs  to  be  computed  numerically 
at  each  time  step  is  the  electric  potential.  The 
current  density  may  be  recast  in  terms  of  the 
effective  electric  field  E '=-Vtp  +  VxB ■  For 
simplicity,  this  expression  is  shown  here  for  a 
magnetic  field  applied  in  the  z-direction  alone: 
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When  the  divergence  of  the  above  expression  is 
set  to  zero  and  it  is  re-written  in  terms  of  the 
electric  potential,  a  Poisson  type  equation  is 
obtained,  which  must  be  solved  at  each 
computational  time  step.  This  equation  is  easily 
derived.  A  full  version  of  this  equation  including 
property  gradients,  is  being  omitted  here  due  to 
space  constraints.  Due  corrections  to  electron- 
beam  current  and  electron  pressure  gradient  will 
be  presented  in  the  final  version  of  the  paper.  At 
interfaces  (solid  as  well  as  liquid,)  where  the 
electrical  conductivity  has  a  gradient,  this 
equation  must  be  written  in  a  weak  form,  since 
the  gradient  of  the  conductivity  is  infinite. 
Virtually,  the  component  normal  to  the  interface, 
of  the  following  equation  must  be  satisfied  in 
order  for  the  solution  to  be  unique: 

V<p  =  VxB 
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Electrical  conductivity  and  the  Hall  parameter 
are  transport  parameters,  that  are  evaluated  using 
appropriate  numerical  expressions,  as  described 
in  Refs.  [2],[  15]. 

Boundary  Conditions: 

Boundary  and  initial  conditions  for  the  equations 
described  must  be  specified  accurately  in  order 
to  prevent  numerical  instabilities  and  non¬ 
physical  solution  behavior.  For  supersonic  and 
super-alfvenic  flows,  information  propagates 
only  in  the  downstream  direction.  The  inflow 
variables  are  then  fixed  and  outflow  variables 
are  obtained  by  extrapolation.  For  subsonic  or 
sub-alfvenic  flows,  characteristic  relations  must 
be  used  in  order  to  propagate  information  either 
into  or  out  of  the  domain. 

No  slip  boundary  conditions  are  applied  to  the 
velocity  along  the  wall  for  viscous  flows.  In  the 
present  abstract,  the  flow  is  considered  to  be 
inviscid,  hence  the  gradient  of  the  tangential 
component  of  velocity  is  set  equal  to  zero  at  the 
walls.  Normal  derivatives  for  species  densities 
are  set  equal  to  zero  on  the  walls.  The  gradient 
of  the  gas  temperature  and  vibrational 
temperature  is  set  equal  to  zero  at  the  walls. 
Electromagnetic  boundary  conditions  may  be 
more  involved,  depending  on  the  choice  of 
material  conditions  at  the  wall.  In  the  case  of  an 
insulating  wall,  the  gradient  of  the  electric 
potential  normal  to  the  wall  is  zero,  and  when 
the  wall  does  conduct,  an  extended  domain  is 
used  to  depict  solid  boundaries.  A  thin 
conducting  wall  option  is  also  available,  wherein 
a  fraction  of  the  incoming  current  is  distributed 
tangentially  along  the  wall,  governed  by  a 
component  of  the  Poisson  equation  for  the 
electric  potential. 

Numerical  Scheme: 

The  governing  equations  are  integrated  using  a 
point-wise  implicit  scheme  described  in 
Munipalli  et  al.  [14],  Upwind  fluxes  based  on 
Roe-type  linearization  are  used  in  a  finite 
volume  format  on  general  unstructured  meshes. 
The  formal  order  of  accuracy  of  this  numerical 
procedure  is  extended  by  using  the  TVD 
approach,  as  described  by  Barth  [3].  Work  is 
underway  to  develop  this  environment  into  a 


higher  order  multi-physics  solver  using  the 
Discontinuous  Galerkin  technique. 

Conducting  walls 

When  the  walls  are  thick,  computations  must  be 
carried  out  within  the  wall  region  for 
electromagnetic  and  thermal  quantities  but  not 
for  fluid  momentum  and  density.  This  involves 
blocking  out  of  these  regions  from  the  flow 
solver.  An  index  array  denotes  the  material  type 
in  each  cell.  Boundary  conditions  are  applied  at 
this  interface  rather  than  the  edges  of  the 
computational  domain,  for  the  flow  solver.  A  far 
field  or  Boundary  Element  Method  (BEM) 
boundary  condition  may  be  applied  at  the  outer 
edge  of  the  solid  walls.  However,  since  solution 
process  must  be  carried  out  internal  to  these 
walls,  it  does  involve  a  sizeable  computational 
expense.  It  is  believed  that  the  BEM  will  reduce 
this  burden  significantly,  and  this  method  is 
described  below.  The  inclusion  of  a  BEM  model 
in  our  MHD  codes  is  presently  ongoing. 

It  is  convenient  to  think  of  the  flow  of  current  in 
MHD  flows  in  terms  of  an  equivalent  static 
electric  circuit,  as  sketched  in  Fig.  [2],  Current 
generated  by  the  flow  of  a  conducting  medium 
in  a  magnetic  field,  flows  in  parallel  in  various 
pathways.  These  include:  (a)  the  conduction  path 
through  the  working  medium-air,  (b)  radiative 
transmission  at  high  power,  charge  separation 
and  cascade  ionization  of  flow,  and,  (c) 
conduction  through  solid  walls  and  other  such 
external  circuits.  As  a  first  cut  estimate,  these 
resistances  may  be  used  to  estimate  the  flow  of 
current  and  the  approximate  Lorentz  force  that 
may  be  generated.  Considering  momentum  and 
energy  balance  across  the  flow  path,  we  may 
estimate  the  pressure  force,  and  from  there, 
control  moments  that  are  produced. 

There  are  three  instances  in  which  the  effects  of 
wall  conduction  become  important  in  MHD: 

(a)  Mesh-resolved  conducting  walls,  where  a 
computational  mesh  is  used  to  resolve  the 
currents  penetrating  the  wall  region.  Fig.  [5] 
shows  a  jet  flow  case  where  current  enters 
the  rectangular  walls  of  a  channel  where  a 
transverse  magnetic  field  is  applied.  Fig.  [3] 
shows  a  case  of  a  channel  wall  with  a  very 
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minute  crack  in  insulation.  Computation 
shows  a  significant  amount  of  current 
leaking  through  this  crack,  causing  a 
pressure  gradient  enhanced  by  about  100 
times  from  the  fully  insulated  case.  These 
results  will  differ  for  different  flow 
geometries  and  MHD  interaction 
parameters.  Hypersonic  flow  cases  will  be 
presented  in  the  final  paper. 

(b)  Infinitely  segmented  walls,  where  an 
assumption  is  used  that  the  current 
tangential  to  the  surface  is  small,  and  that 
the  electric  potential  gradient  in  the 
tangential  direction  may  be  specified 
analytically.  Figure  [4]  shows  a  case  of  an 
MHD  accelerator  with  segmented 
electrodes,  where  prohibitively  large  mesh 
sizes  are  needed  to  resolve  individual 
electrodes.  Analytical  BCs  provide  a  simple 
means  to  compute  the  flow. 

(c)  Thin  conducting  walls,  where  a  continuity 
equation  for  current  (in  terms  of  the  electric 
potential)  is  written  in  a  wall-tangential 
direction.  This  model  has  been  developed 
for  incompressible  flow  applications  at 
HyPerComp  (Refs.  [16],  [17]),  and  is 
presently  being  added  to  the  hypersonic 
flow  solver. 

Coupled  Boundary  Element  Method 
The  Boundary  Element  Method  (BEM)  is  an 
attractive  choice  for  some  of  the  coupled  solid- 
fluid  problems  encountered  in  MHD.  Simply 
put,  this  method  reduces  a  volume  problem  to  a 
surface  problem,  thereby  reducing  its 
dimensionality  by  one,  and  thus  avoiding  the 
need  to  generate  a  volume  mesh  in  semi-infinite 
or  completely  enclosed  regions.  Examples  could 
be  the  air  enclosing  the  outer  walls  of  MHD 
chambers,  metallic  or  insulating  penetrations  in 
liquid  metal  flow,  the  far  field  in  a  two  fluid 
flow  and  so  forth.  The  importance  of  BEM  has 
been  recognized  by  various  authors.  Masse  et  al. 
[10]  and  related  codes  at  http://www. 
simulog.fr/eis/fluxl  ,htm.  Tezer-Sezgin  et  al. 
[21]  and  Guermond  et  al.  [5]  present  a  flavor  of 
the  potential  of  this  method.  A  possible 
application  of  this  method  is  presented  here,  that 
can  be  used  for  both  MHD  as  well  as  heat 
transfer  applications  in  much  the  same  manner. 


The  most  significant  use  of  BEM  is  seen  in 
computational  regions  where  there  is  no  flow 
(Ref.[21]  uses  the  BEM  for  the  flow  itself,  but 
with  certain  simplifying  assumptions.)  In  a 
steady  state  situation,  both  heat  conduction,  the 
diffusion  of  magnetic  field  as  well  as  the 
electrical  potential,  all  obey  Laplace’s  equation, 
which  can  be  quickly  written  in  a  discrete 
boundary  dependent  manner:  For  an  arbitrary 
quantity  u,  we  have: 


The  starred  quantity  u*  represents  a  Green’s 
function  for  Laplace’s  equation,  that  is 
analytically  known.  Boundary  conditions  may 
prescribe  u  or  the  normal  derivative  of  u  or  a 
linear  combination  of  both. 


Results 


In  this  paper,  we  present  results  obtained  by 
solving  for  the  electromagnetic  fields  within  the 
solid  walls,  in  this  case,  using  the  electric 
potential  formulation.  Some  frequently  used 
non-dimensional  quantities  in  this  type  of  study 
are: 


<j  t  , 

Wall  conductance:  c  —  — — -  -  the  ratio  of  the 
a0L 

wall  conductivity  to  the  mean  flow  conductivity, 
multiplied  by  ratio  of  the  wall  thickness  to  the 
characteristic  flow  length. 


Hartmann  number: 


Ha  =  BL 


where  B  is 


the  mean  magnetic  field  applied. 

Interaction  parameter:  N  = - ,  which  scales 

pV 

the  effect  of  the  MHD  interaction  with  flow 
momentum. 


Using  the  estimates  for  the  resistances  of  the 
various  parts  of  the  flow  described  in  Figure  [2], 
Tillack  [21]  derived  a  method  to  evaluate 
pressure  drop  across  an  MHD  channel  with  fully 
developed  flow  and  conducting  walls.  Figure  [] 
uses  this  analysis  to  estimate  pressure  gradient  in 
a  fully  developed  flow  at  various  Hartmann 
numbers  for  various  values  of  wall  conductance, 
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at  a  Reynolds  number  of  106.  Figure  []  shows  a 
sample  computation  along  these  lines  of  a 
developing  flow  with  conducting  walls.  It  is 
rather  apparent  that  the  effect  of  wall 
conductivity  on  flow  pressure  gradient  is 
sizeable.  It  must  be  noted  that  these  estimates 
are  based  on  incompressible  flow  assumption. 
However,  the  situation  is  not  vastly  different  for 
compressible  flows. 

We  present  now,  a  different  look  at  this  situation 
by  considering  external  flow  past  a  cylinder.  In 
ideal  fluid  dynamics,  viscous  flow  past  a 
cylinder  produces  the  Karman  vortex  street 
downstream,  and  an  oscillatory  pressure  field  in 
the  aft  side  of  the  cylinder.  When  mean  flow 
quantities  are  considered,  statistical  averages  of 
the  recirculation  length  behind  the  cylinder,  and 
other  such  global  effects  may  be  discerned. 
Figure  [7]  shows  the  mean  pressure  field  in  the 
evolution  of  the  flow  past  a  cylinder  with  and 
without  MHD  and  conducting  walls.  A  cursory 
glance  reveals  that  the  pressure  drop,  thereby 
drag,  across  the  cylinder  is  reduced  with  MHD, 
due  to  the  relative  uniformity  of  pressure  in  the 
close  vicinity  of  the  object.  An  additional  feature 
of  interest,  is  that  this  effect  is  more  prevalent  in 
the  insulating  wall  case  rather  than  the 
conducting  wall  case,  because  of  the  manner  in 
which  current  is  distributed  in  the  flow  field  in 
both  cases.  The  current  lines  are  shown  in 
Figure  [6]. 

We  present  two  cases  of  supersonic  flow  in  a 
channel  with  conducting  walls.  In  the  first  case, 
the  walls  are  flat,  and  the  flow  is  smooth  inside 
the  channel.  The  MHD  effect  is  weak,  the 
interaction  parameter  being  rather  small.  The 
main  reason  for  this  is  that  the  conductivity  is 
rather  small,  and  the  three  dimensional  nature  of 
MHD  has  not  been  used  to  advantage.  Indeed, 
the  magnetic  field  is  applied  in  the  direction 
normal  to  the  plane,  and  currents  there  fore  must 
flow  within  the  plane  of  the  flow.  First,  we  study 
the  amount  of  current  that  is  generated  within 
the  flow  for  conducting  as  well  as  insulating 
walls.  Figure  [8a]  shows  the  case  when  the  flow 
conductivity  is  about  100  mho/m,  while  the  wall 
conductivity  is  l.e-6  mho/m.  Figure  [8b]  shows 
the  roles  reversed,  when  the  wall  conductivity  is 
l.e3  mho/m,  and  the  flow  conductivity  is  about 


10  mho/m.  The  difference  in  the  magnitude  of 
the  current  that  is  flowing  through  the  channel  is 
phenomenal,  being  a  few  orders  of  magnitude 
higher  for  the  case  with  conducting  walls,  on 
average.  If  the  interaction  parameter  is  increased 
for  this  calculation,  no  stable  solution  was 
observed  numerically.  This  is  attributed  to  the 
electron  heating  model  currently  used  in  the 
code,  that  seems  to  cause  rather  violent  shock 
disturbances  caused  by  Joule  heating.  A  more 
gradual  relaxation  model  between  electron 
energy  and  flow  energy  is  being  incorporated. 
Current  lines  are  sketched  in  Figure  [9]  for  both 
these  cases. 

Figure  [10]  shows  the  flow  in  a  supersonic 
compression  channel  with  conducting  walls  at  a 
Mach  number  of  3.  Again,  the  interaction 
parameter  was  held  deliberately  low,  just  to  get  a 
feel  for  the  flow  currents  that  are  produced.  We 
same  feature  was  again  seen,  namely,  that  the 
conducting  wall  case  has  about  10  times  higher 
current  in  the  flow  on  an  average  than  the 
insulating  wall  case.  Figure  [11]  shows  the 
shock  pattern  in  this  channel. 

We  are  currently  attempting  to  integrate  these 
models  with  other  advanced  physical  models  for 
ionizing  air  for  realistic  hypersonic  flight 
conditions.  Figure  [12]  shows  sample  plots  of 
current  lines  from  our  recent  study  in  hypersonic 
inlets  with  insulating  walls.  The  same  case  will 
be  applied  to  conducting  walls  in  the  near  future. 

Conclusion 

This  paper  presents  the  results  of  our  ongoing 
investigations  in  the  development  of  realistic 
physical  and  numerical  models  in  plasma  and 
MHD  flows.  We  have  seen  the  effects  of  wall 
conductivity  to  be  crucial  to  pressure  variations 
in  both  internal  and  external  flows,  and  have  an 
impact  on  drag,  mass  capture  and  other  such 
parameters  affecting  hypersonic  flow  systems. 
In  the  future,  we  intend  to  demonstrate  these 
effects  on  a  full  scale  hypersonic  NASP  type 
inlet  configuration. 
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Insulator 


Figure  -2  :  A  circuit  diagram  for  current  flow  in  MHD 
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Figure  4  :  Current  lines  in  the  vicinity  of  a  pair  of  powered  electrodes  (above,)  and 
Potential  distribution  along  an  MHD  accelerator  channel  powered  by  segmented  electrodes 
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-  dp/dx 


Figure  5:  Closed  channels  with  conducting  walls  at  Re  =  106,  showing  pressure  drop  and  current  vectors 
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Figure  6:  Cylinder  flow  with  and  without  conducting  walls 
Ha  =  1 0,  c  =  20  (left),  c  =  0.002  (right) 
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Figure  7:  Clockwise  from  left:  (a) 
pressure  contours  without  MHD,  range  - 
0.8  to  0.4,  (b)  with  conducting  walls,  c  = 
20,  range  -0.4  to  0.4,  (c)  with  insulating 
walls,  c  =  0.002,  range  -0.4  to  0.4 
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Figure  8:  Supersonic  channel  flow  with  a  fringing  B-field,  insulating  walls  (left)  conducting  walls  (right) 


Figure  9:  Current  lines  for  the  above  cases 
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Figure  10:  Current  magnitude  for  a  channel  flow  in  an  inlet  type  configuration  with 
insulating  walls  (left)  and  conducting  walls  (right) 
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